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AN  OVERVIEW  OF 


RECURSIVE  LEAST  SQUARES  ESTIMATION 
AND  LATTICE  FILTERS  t 


John  M.  Turner 


1.  INTRODUCTION 


Tbe  lattice  liter  streets  re  is  an  alternative  means  of  realising  a  digital  Biter  transfer  func¬ 
tions.  Although  Ike  lattice  liter  structure  (also  called  the  ladder  structure)  does  not  have  the 
m.awwm  umber  of  multipliers  and  adders  for  a  transfer  function  realisation,  it  does  have  several 
advantageous  properties.  These  include  cascading  of  identical  sections,  coefficients  with  magni¬ 
tudes  less  than  one,  stability  test  by  inspection,  and  good  numerical  roundoff  characteristics. 
Moreover,  tbe  lattice  liter  structure  is  particularly  suited  for  adaptive  lltering  since  the  recursive 
solution  of  least  squares  estimators  naturally  produces  n  lattice  liter  structure.  Also  tbe  lattice 
liter  structure  octbogouafises  the  input  signal  on  a  stage  by  stage  basin  This  leads  to  very  fast 
convergence  and  tracking  capabilities  of  the  lattice  structure.  Although  many  alternative  tech¬ 
niques  have  been  developed  to  estimate  the  relectios  coefficients  that  parameterize  the  lattice 
structure,  the  recursive  least  squares  method  updates  the  least  squares  estimate  upon  the  observa¬ 


tion  of  each  data  sample.  This  procedure  leads  to  an  optimal  estimate  and  requires  only  slightly 
greater  computational  burden  than  alternative  techniques. 

~^Thii  overview  presents  the  derivation  of  the  recursive  least  squares  lattice  liter  using  n 
recursive  extension  of  the  standard  block  data  Levinson  least  squares  solution.  The  linear  predie- 
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tor  liter  presented  here  has  been  widely  used  for  synthesis  of  speech  waveforms  (LPC),  deconvo¬ 
lution  of  seismic  data,  high  resolution  spectral  estimation,  adaptive  line  enhancement,  adaptive 
noise  cancelling  and  adaptive  antenna  array  processing.  The  ideas  embodied  in  this  estimation 
technique  are  derived  from  the  work  of  many  individuals  since  1970.  The  approach  presented 
here  follows  that  of  [Lee,  1980]. 

Adaptive  estimation  techniques  modify  the  estimation  filter  parameters  according  to  the 
newly  observed  data  sample.  For  every  new  data  sample,  recursive  estimation  using  the  lattice 
filter  generates  new  reflection  coefficients  and  prediction  errors  for  every  filter  order.  Changing 
every  filter  coefficient  for  each  new  data  sample  is  important  for  applications  where  fast  conver¬ 
gence  or  tracking  of  quickly  time  varying  signals  is  required.  However,  for  applications  where  the 
dynamics  are  slow,  only  the  results  after  observing  the  signal  for  a  certain  time  period  are  impor¬ 
tant.  The  recursive  algorithms  described  here  can  also  be  used  to  accumulate  signal  properties 
over  a  particular  time  period.  The  procedure  for  converting  between  lattice  filter  coefficients  and 
the  more  common  equivalent  tapped  delay  line  filter  coefficients  is  given  in  Section  2. 

The  mathematics  of  recursive  least  squares  estimation  requires  the  updating  of  variables 
with  time  and  order  subscripts.  The  algorithms  for  this  are  often  complicated.  Therefore,  an 
intuitive  introduction  in  Section  3  presents  the  nature  of  the  lattice  filter  structure,  the  stage  by 
stage  orthogonalising  property,  and  analogies  with  physical  phenomena.  Section  4  briefly  presents 
approximation  techniques  for  determining  the  reflection  coefficients  from  observed  data.  The 
advantage  of  the  lattice  filter  structure  is  that  time  recursive  exact  least  squares  solutions  to  esti¬ 
mation  problems  can  be  efficiently  computed.  The  development  of  the  Recursive  Least  Squares 
Lattice  estimation  algorithm  is  presented  in  Section  5  and  8.  A  square  root  normalised  least 
squares  lattice  algorithm  that  has  better  numerical  properties  is  presented  in  Section  7. 

The  computational  complexity  of  these  algorithms  is  discussed  in  Section  8.  An  efficient 
means  of  implementing  the  recursive  least  squares  algorithm  using  rotational  arithmetic  is 
presented.  This  rotational  arithmetic,  called  CORDIC  arithmetic,  is  not  new,  having  been  used 
for  calculating  trigonometric  Auctions  in  hud  held  calculators.  The  design  of  u  integrated  cir- 
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cuit  chip  to  implemeBt  the  least  squares  lattice  algorithm  using  CORDIC  arithmetic  is  mentioned. 

To  demonstrate  the  power  of  this  adaptive  estimator,  simulation  examples  of  convergence 
and  tracking,  examples  using  real  speech  data  and  electrophysiologieal  data  and  adaptive  equal¬ 
iser  examples  are  presented  in  Section  9.  Since  recursive  least  squares  estimation  and  lattice  filter 
structures  have  been  a  very  active  area  of  research,  Section  10  refers  to  related  ideas. 

Since  the  equations  developed  here  are  recursive  in  order  and  time,  the  following  notation  is 
used.  A  variable  *(<)  is  a  general  time  sampled  data  value  while  »T  is  the  specific  data  sample  T 
samples  after  the  beginning  of  the  recursion  (relative  time  T).  Bold  capital  letter  variables 
represent  matrices  or  vectors.  When  two  subscripts  are  used,  the  first  is  the  order  and  the  second 
is  the  time  parameter,  ie.  Aj  t  is  the  vector  of  i-tk  order  predictor  coefficients  determined  from 
data  up  to  the  specific  time  T. 
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I.  GENERAL  LATTICE  DIGITAL  FILTER  STRUCTURE 


The  general  lattice  digital  titer  is  a  mease  of  realizing  a  digital  titer  transfer  function.  The 
lattice  s  tracts  re  itself  is  motivated  by  similar  analog  titer  structures  that  have  good  properties. 
In  this  section,  the  digital  lattice  titer  is  introduced  and  related  to  the  direct  form  tapped  delay 
line  digital  titer. 

Since  analog  lattice  and  ladder  titers  have  desirable  characteristics,  digital  titers  with  simi¬ 
lar  structures  were  investigated.  For  example,  third  order  LC  Butterworth  titer  shown  in  Fig.  la 
is  a  simple  analog  lattice-type  structure.  It  is  noted  for  the  relative  insensitivity  of  its  frequency 
response  to  slight  perturbations  of  the  circuit  element  values  around  their  nominal  values.  This 
structure  can  be  transformed  into  the  general  lattice  network  shown  in  Fig.  lb.  The  latter  is  the 
lattice  structure  of  interest  in  this  chapter.  The  analog  lattice  structure  consists  of  a  cascade  of 
identical  stages,  each  stage  with  a  pair  of  input  and  output  terminals.  By  developing  a  digital 
liter  conlga ration  that  is  similar  to  the  analog  lattice  structure,  the  digital  liter  inherits  many  of 
the  same  properties.  Since  the  structure  of  a  digital  liter  realisation  inluences  its  sensitivity  to 
Inite  word  length  arithmetic,  the  digital  lattice  liter  has  good  numerical  properties. 

The  digital  lattice  liter  realisations  consist  of  cascaded  stages  with  two  input  and  two  out¬ 
put  ports,  as  in  the  analog  structures.  Possible  digital  lattice  configurations  for  realising  a  general 
digital  transfer  function  include:  an  asymmetric  multiplier  form  (Fig.  2)  and  a  symmetric  two 
multiplier  form  (Fig.  3).  For  the  asymmetric  multiplier  lattice,  the  structure  inside  each  stage 
realms  a  single  pole  and  zero  equivalent  transfer  function.  The  algorithm  for  determining  the 
asymmetric  multiplier  structure  of  Fig.  2  can  be  found  in  (Kfitra  et  al,  1077].  This  form  degen¬ 
erates  to  a  tapped  detap  Has  for  either  all  pole  or  all  aero  transfer  functions  and  thus  is  not  of 
interest  here.  The  symmetric  two  multiplier  form  does  not  degenerate  but  it  requires  more  multi- 
plots  than  an  equivalent  tapped  delay  lae  liter.  This  lattice  liter  can  be  modified  to  a  one  mul¬ 
tipier  form  so  as  to  have  the  minimum  number  of  multi  pliers,  but  this  requires  extra  adders  (Fig. 
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A  cascade  of  lattice  sections,  forming  a  lattice  filter  can  implement  a  digital  transfer  fane* 
tion  in  n  way  that  has  advantages  over  the  standard  direct  form,  parallel,  or  standard  cascade 
form  realisations.  The  cascaded  structure  in  the  lattice  filter  propagates  a  forward  signal  f  j(  t ) 
and  a  backward  signal  ly(<)  at  time  t  and  section  number  j.  The  fundamental  equation  describ> 
ing  the  lattice  filter  structure  is  (1),  see  Figure  4. 

/j+i(0-//(0-*/+i*iM)  (i) 

*/+ 1(0  -  M*-1)  -  */+ 1  /AO 

The  multiplieis  in  the  crossover  portion  of  the  lattice,  k  are  known  as  reflection  coefficients  or 
partial  correlation  (PAR COR)  coefficients. 

The  implementation  of  digital  filter  transfer  functions  in  lattice  form  have  been  examined 
(Gray  and  Market,  1973,1075].  State  space  canonical  forms  were  also  established  [Morf,  1974, 
Morf  et  al,  1977,  Lee,  1980].  ALGORITHM  1  determines  the  reflection  coefficients  £,  and  tap 
coefficients  v,  for  the  lattice  filter  of  Fig.  3  that  is  equivalent  to  a  (stable  nonreducible)  direct 
form  transfer  function  with  numerator  coefficients  if  and  denominator  coefficients  §f  (from  (Gray 
and  Market,  1973)).  The  one  multiplier  form  in  Fig  3b  uses  coefficients  v’.  Although  the  lattice 
coefficients  and  the  direct  form  coefficients  are  related  in  a  nonlinear  manner,  this  algorithm  is 
invertible  so  that  a  lattice  structure  can  be  converted  uniquely  to  a  direct  form  filter  and  vice 
versa  (when  aO  the  roots  are  inside  the  unit  circle). 


ALGORITHM  1: 
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General  Transfer  Function  to  Lattice  Filter 


E  */  ** 

HA*)  -  nr -  ;  ••  - 1 

e«; 

>-« 


*r  “i 

For  »  —  P  to  I 

»<  —  */ 

*  *  —  *i  /  *< 

*i-i  *■*<(!+*<) 

For  /  —  1  to  i-1 


'  1  - 
*r  -  *i  -  %  •/-/ 


continue 
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continue 
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END 

While  Fig.  3  and  ALGORITHM  1  describe  the  lattice  filter  for  a  general  transfer  function 
poles  sad  seros,  the  remainder  of  thin  chapter  discusses  all  pole  transfer  functions  and  their 


with  poles  and  seros,  the  remainder  of  this  chapter  discusses  all  pole  transfer  functions  and  their 
inverses,  al  sera  transfer  tactions.  For  an  aO  pole  transfer  function  (i£  ■  1,  if  ■  0,  ;>1),  the 


lattice  filter  is  callsd  the  feedback  lattice  liter  shown  in  Fig.  4.  The  inverse  of  the  feedback  lat¬ 
tice  can  be  determined  by  applying  Mason’s  rule  to  Fig.  4.  This  finite  impulse  response  filter,  an 
al  sero  transfer  function,  is  the  feedforward  lattice  (Fig.5).  Thus  a  feedforward  lattice  and  a 
feedback  lattice  with  the  same  coefficients  perform  inverse  operations  on  the  input  signal.  If  a 
signal  is  spphsd  to  a  feedforward  lattice  filter  and  the  result  is  applied  to  a  feedback  lattice  filter, 
the  original  signal  is  retimed.  According  to  Mason’s  rule,  the  reflection  coefficients  parameterize 
hath  the  feedback  and  feedforward  lattice  with  the  appropriate  change  in  signal  low.  ALGO¬ 
RITHM  2  gives  the  procedure  for  converting  horn  reflection  coefficients  to  tapped  delay  line 


the  signal  flow 


whether  the  all  pole  or  all  sero  transfer  function  is 
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ALGORITHM  2: 


Lattice  Coefficients  to  Tapped  Delay  Line  Coefficients 

s* 

For  *  —  2  to  P 

tf  mm  —  kj 

For  j  mm  1  to  »-l 

.f  _  4-1  j,  .i-1 

•y  *  Sy  tj  Sj-y 

continue 

continue 

END 

The  coefficient  sensitivity  of  the  lattice  implementation  of  a  general  digital  transfer  function 
has  not  been  studied  as  thoroughly  as  other  common  filter  structures.  Scaling  conventions 
[Market  and  Gray,  1975(2)]  and  roundoff  noise  characteristics  [Markel  and  Gray,  1075(1)]  for  finite 
wordlength  arithmetic  were  developed  for  various  lattice  configurations.  The  one  multiplier  lat¬ 
tice  has  roundoff  none  characteristics  that  are  always  better  than  the  two  multiplier  lattice.  Both 
lattice  filters  are  always  better  than  the  direct  form  realisation.  Particularly  when  the  width  of 
the  filter  pass  band  becomes  small,  the  The  lattice  liter  structure  has  better  roundoff  noise 
characteristics  than  other  filter  realisations,  particularly  when  the  width  of  the  filter  pass  band 
becomes  small.  A.  normalised  lattice  filter  (requiring  more  multiplications)  was  developed  that 
performs  better  than  the  other  lattice  structures  or  parallel  form  realisations. 

The  implementation  of  a  transfer  function  requires  quantised  coefficients  that  can  effect  the 
stability  of  a  filter  and  its  inverse.  The  sensitivity  of  the  roots  of  the  transfer  function  to  pertur¬ 
bations  of  the  lattice  filter  and  tapped  delay  line  filter  coefficients  has  been  investigated  [Chu  and 
Messerschmitt,  1980,1983).  The  effect  of  varying  the  tapped  delay  line  filter  coefficients  was  the 
same  lot  each  coefficient.  When  the  roots  are  dose  to  the  unit  circle,  quantisation  of  the  tapped 
delay  fine  coefficients  tends  to  move  the  roots  perpendicular  to  the  unit  circle.  For  an  all  pole 
transfer  fraction,  this  quantisation  can  cause  the  polm  to  move  outside  the  unit  circle  and  the 
transfer  function  to  become  unstable.  For  lattice  liters,  the  effect  on  the  root  location  of  varying 
the  tow  order  coefficients  is  :  ueh  greater  than  for  the  higher  order  coefficients.  For  relection 
eeeffideote,  the  rwr  tm.  to  more  tangentially  to  the  unit  circle,  thus  changing  the  center 
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frequency  rather  than  the  bandwidth  of  the  roots.  Low  order  reflection  coefficients,  particularly 
those  with  magnitudes  near  unity  need  to  be  more  accurately  quantized  than  the  higher  order 
reflection  coefficients.  No  simple  rule  of  thumb  exists  for  tapped  delay  line  filters.  The  effects  of 
the  lattice  coefficient  quantisation  have  been  studied  most  extensively  for  speech  modeling  appli¬ 
cations.  For  typical  prediction  filters  used  in  speech  processing,  substantially  coarser  quantization 
of  the  reflection  coefficients  than  of  the  tapped  delay  line  filter  coefficients  is  possible  while  still 
maintaining  the  subjectively  perceived  spectral  response. 


C  ^  C 


Fig.  la  Third  Order  LC  analog  Ladder 


Fig.  lb  General  Analog  Impedance  Lattice 
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3.  PROPERTIES  OF  THE  LATTICE  STRUCTURE 


The  lattice  filter  has  a  more  complex  structure  and  requires  more  numerical  operations  to 
implement  a  transfer  function  than  the  straight  forward  direct  form  realization.  However  this 
increased  complexity  is  offset  by  several  advantageous  properties  of  the  lattice  structure  including; 
a  stability  test  by  inspection,  a  stage  by  stage  orthogonalization  of  the  input  signal,  and  a  physi¬ 
cal  interpretation  as  wave  propagation  in  a  stratified  medium.  The  lattice  filter  structure  natur¬ 
ally  evolves  from  a  prediction  filter  where  orthogonality  conditions  are  applied.  The  properties  of 
the  lattice  estimation  filter  are  presented  in  thia  section.  An  acoustical  tube  model  of  the  human 
vocal  tract  is  interpreted  as  a  lattice  filter  in  Section  3.2  to  lend  a  physical  interpretation  to  the 
reflection  coefficients. 


3.1  Orthogonarilxlng  Properties 


In  early  investigations  of  the  lattice  structure,  a  connection  with  orthogonal  polynomials  was 
noted  [Itakura  and  Saito,  1968,  Burg,  1975,  Makhoul,  1975,  Market  and  Gray,  1976].  Lattice 
form  realisations  are  obtained  by  an  orthogonalization  of  the  state-space  transfer  function  using 
Szegp  orthogonal  polynomials.  The  theory  of  Ssego  polynomials  and  their  applications  in  system 
theory  (stability  testing)  and  in  stochastic  problems  (prediction  theory  and  spectral  analysis)  has 
been  discussed  in  (Grenander  and  Szego],  The  Schur  stability  test  uses  the  properties  of  orthogo¬ 
nal  polynomials  to  determine  whether  the  poles  of  a  transfer  function  are  inside  the  unit  circle 
and  hence  a  stable  transfer  function.  The  test  is  performed  by  using  ALGORITHM  1  to  compute 
the  lattice  filter  then  checking  that  the  magnitudes  of  all  the  reflection  coefficients  {£,}  are  less 
than  one. 


For  problems  in  estimation,  the  minimum  mean  squared  error  estimation  equations  can  be 
transformed  into  a  stage  by  stage  optimisation  by  forming  a  recursion  on  the  optimum  filter 
order.  The  parameter  in  the  recursion  can  be  estimated  stage  by  stage  since  it  depends  on  quanti¬ 
ties  that  are  orthogonal  between  stages.  This  orthogonalization  property  is  now  developed  for  the 


prediction  filter  (an  all  tero  transfer  function). 

For  a  predictor  of  p-th  order,  the  data  sample  at  time  f,  x(<),  is  approximated  as  a  linear 
combination  of  the  previous  p  samples,  x(<-l),...,x(i-p).  The  forward  prediction  error  /,(*) 
must  be  orthogonal  to  the  previous  data  samples  to  attain  the  minimum  mean  squared  error 
value.  This  determines  the  weighting  factors,  {a,}  on  the  previous  data. 

f,(t)  —  x(l)  +  s,x(<-I)  +  •  +  s,x(f-p)  (2) 

x(t-j))  -  0,  1  <j<p 

The  operation  E(  )  represents  the  statistical  expectation.  This  prediction  error  is  also  called  a 
prediction  residual  or  an  innovation  when  the  coefficients  are  chosen  to  attain  the  minimum  mean 
squared  error.  A  backward  prediction  error,  can  similarly  be  defined  to  predict  x(t-p-l) 

from  the  same  samples,  x(f-l),...,s(f-p). 

M*"1)  “  *(<-*-!)  +  «i *0-p)  +  •  *  +  crz(t- 1)  (3) 

*(!-/))  -  0,  I <j<p 

Here  the  {c7}  are  chosen  to  satisfy  this  condition.  Notice  that  both  prediction  errors  satisfy  the 
same  orthogonality  conditions. 

Increasing  the  prediction  order  to  y-f  1,  V  i(f )  represents  the  component  of  x(f)  that  is  not 

predictable  from  z(f-l) . x(l-y),x(f-p-l).  The  p-th  prediction  error  uses  information  up  to 

x(f-p),  so  now  the  information  about  x(f )  that  can  be  predicted  from  x(f-p-l)  must  be  included. 
However  much  of  this  information  is  already  contained  in  z(f-l),...,x(f-p).  The  backward  pred¬ 
iction  error  represents  the  new  information  in  the  sample  z(f-p-l).  The  plausible  recur¬ 

sion  for  fp+  i(t )  is  (4)  where  the  scalar  x  is  determined  so  that  i(t)  satisfies  the  new  ortho¬ 
gonality  conditions. 

/,♦!<«) -/,(*)-*/♦»  MM  (4) 

£(/,♦»(*)  *<M)  -  0,  l£;£y+i 

The  only  constraint  not  immediately  satisfied  Involves  s(f-p-l),  and  is  given  by  (5).  By  substi- 


toting  (3)  and  (4)  in  (S),  the  optimal  1  ia  detennined  (6). 


(5) 

0  -  E[J9(t)  x(t-p- 1))  -  4/+1  E(4,(t-1)  *(l-y-l)) 

-E(f,(t)k,(t-l))-k}+l  £(*,*(1-1)) 

*/+  i  (e) 

Similarly,  the  recnnion  for  the  backward  predictor  is  obtained  (7)  and  the  optimal  x  is  deter* 

mined  (8). 

*,+  i(0  —  *,(*-!)-  *>+i  /»(*)  (7) 

*(0)-o 

0  -E(f,(t)i,(t-l))-k>+ iE(f,*(t)) 

*,*>!  -  E(f,0)  /  BUM  («) 

Extending  the  prediction  liter  to  the  next  higher  order,  y+  2  requires  the  calculation  of  the 
new  prediction  errors,  /,♦!  and  *,♦!  from  (4)  and  (7).  Thus  a  prediction  liter  can  be  constructed 
solely  using  the  lattice  structure  by  successively  increasing  the  liter  order.  This  is  the  stage  by 
stage  orthogonaliiation  property  of  the  lattice  structure  where  each  relection  coefficient  is  deter¬ 
mined  separately.  This  stage  by  stage  computation  of  prediction  coefficients  does  not  hold  for  the 
tapped  delay  line  liter  (2).  The  coefficients  {*;}  are  interdependent  and  they  all  change  when 
the  liter  order  increases. 

Further  insight  into  properties  of  the  prediction  errors  is  provided  in  (Makhoul,  1978(2)]. 
The  backward  prediction  error  results  from  a  Gram-Schmidt  type  orthogoaalisatioo  of  delayed 
versions  of  the  signal  This  property  of  orthogonal  variables  makes  the  lattice  structure  advanta¬ 
geous  for  adaptive  lltering.  Also  the  decrease  ia  signal  energy  after  each  prediction  stage  is  easily 
determined.  This  feature  eaa  be  used  to  scale  the  prediction  errors  to  maintain  good  numerical 
properties.  The  most  important  properties  are  summarised  here. 


(») 


E(/,(t)  */(<-!))  "  |q  l£j<p 


£(/,*(<))  -  £(/,(*)*(<))-*/ 

£(/,(<)  />(*))  -  *}  i<J<p 


(10) 


E(b*it-l))  -  F(i,(l-l)z(l-p-l))  -  *i  (11) 

E(»,((-l)  M‘-»)  -  (o'  [</<P 

*1+ 1  -  «/  ( 1  -  V  V  )  (12) 

«*♦»  -»/(!-*/  *») 

When  the  signal  z()  is  stationary  with  known  autocorrelation  function,  the  forward  and 
backward  prediction  error  energies  at  each  stage  are  identical  (o/  —  <rf).  Then  the  two  reflection 
coefficients  are  equal  and  the  symmetric  two  multiplier  lattice  structure  computes  these  prediction 
error  recursions.  When  the  signal  to  be  modeled  is  assumed  to  be  stationary,  a  single  reflection 
coefficient,  k  is  determined  by  combining  sample  data  estimates  of  t*  and  1*.  This  lattice  filter 
with  constant  coefficients  is  the  feedforward  lattice  (1)  of  Section  2.  For  nonstationary  signals, 
adaptive  estimates  are  generated  by  making  the  reflection  coefficients  time  varying. 

The  reflection  coefficients  are  closely  related  to  partial  correlation  factors  which  have  several 
interesting  statistical  properties.  The  correlation  between  t(t)  and  z(l-p-l),  after  their  mutual 
linear  dependence  on  the  intervening  samples  (z(t-l) , . . . ,  z(t-y)}  has  been  removed  is 
£(/,(!)  i,(i)).  This  relation  arises  from  the  orthogonalixiag  nature  of  the  lattice.  When  this 
correlation  is  normalised  by  the  variance  of  /,  and  bf,  it  is  known  as  the  p-th  order  partial 
correlation.  The  autocorrelation  function  of  a  stationary  unit  variance  discrete  time  process  can 
be  uniquely  characterized  by  a  sequence  of  reflection  coefficients,  having  values  less  than  or  equal 
to  one  [Barndorff- Nielsen  and  Schou,  Ramsey],  For  any  p-th  order  AR  process,  the  partial  corre¬ 
lation  of  higher  order,  lag  y  -f  i,  (i  *■  1,  2,  . .)  is  aero.  For  a  stationary  AR  process,  the  sample 
estimates  of  the  partial  correlations  are  asymptotically  Ganssian  and  independent  (see  (Mnrthy 
and  Narastmhamj  for  more  statistical  properties). 

b  appications  sack  as  noise  caacelliag  or  equalisation,  the  ortkogoaalising  properties  of  the 
lattice  are  of  primary  interest  to  obtain  fast  tracking  or  convergence,  see  Section  0.  The  back* 


ward  prediction  errors,  i,(f)  are  extensively  used  since  they  are  a  Gram-Schmidt  orthogonality 
tion  of  delayed  versions  of  the  input  time  series. 

U  Physical  Interpretation 

The  lattice  structure  and  the  reflection  coefficients  have  a  physical  interpretation  that  for 
particular  classes  of  signals  lends  understanding  to  the  properties  of  the  lattice  structure.  Model¬ 
ing  of  wave  propagation  in  a  stratified  medium  leads  to  a  cascade  of  lattice  filters.  Thu  model 
has  been  applied  in  seismic  signal  processing  by  [Treitel  and  Robinson,  Burg,  1067]  and  others. 
The  physical  properties  of  scattering  medium  leads  to  inversion  methods  based  on  cascaded 
reflection  elements,  eg.  the  characterisation  of  (electrical)  transmission  lines  [Gopinath  and  Son- 
dhi,  1971]  or  the  human  vocal  tract  [Gopinath  and  Sondhi,  1970).  Similarly  in  the  fields  of  acous¬ 
tics  and  speech  processing,  an  acoustic  transmission  line  with  step  changes  in  impedance  leads  to 
a  lattice  cascade  structure.  The  human  vocal  tract  has  been  modeled  as  a  cascade  of  acoustic 
tube  sections  with  different  impedances.  This  relationship  between  a  physiological  system  and  the 
lattice  structure  gives  a  physical  meaning  to  the  reflection  coefficients  and  led  to  the  development 
of  speech  synthesis  systems  using  the  lattice  structure.  The  remainder  of  this  section  develops  an 
acoustical  tube  model  of  the  vocal  tract  into  a  lattice  filter  (see  [Flanagan,  Market  and  Gray, 
1976,  Rabiner  and  Schafer,  1978]). 

A  lossless  acoustical  tube  transmission  line  composed  of  cascaded  cylinders  of  differing  diam¬ 
eter  but  equal  length  was  developed  as  a  model  of  the  vocal  tract  in  {KeQy  and  Lochbaum).  This 
vocal  tract  model  was  studied  to  obtain  a  better  understanding  of  the  speech  production  mechan¬ 
ism  and  to  synthesise  speech  by  computer.  Speech  sounds  result  from  pressure  waves  resonating 
in  the  vocal  tract  (acoustic  tube).  The  significance  of  the  model  is  that  the  cascaded  cylinders 
become  cascaded  lattice  stages.  The  cross  sectional  areas  of  adjacent  cylinders  specify  reflected 
and  transmitted  acoustic  wave  components  which  translate  into  the  lattice  reflection  coefficients. 

Sound  waves  that  propagate  in  a  cylindrical  section  obey  the  conservation  of  momentum 
and  mam  equations  (assuming  standard  conditions,  see  [Rabiner  and  Schafer,  1978j).  Since  the 


cross-sectional  ire*  of  the  n-th  tube  is  constant,  combining  the  conservation  laws  yields  a  one 
dimensional  wave  equation.  To  satisfy  this  equation,  the  steady  state  volume  velocity  «<(*,!)  and 
the  pressure  wave  y,(*,l)  are  composed  of  waves  traveling  in  the  forward,  »+  and  backward,  u~ 
direction. 


dp  p  du 
di  “  “  A  dt 

dx  pea  dt 

%(».*)  —  ui+(t-s/c)-ui-{t+x/c)  (13) 

Piix.t)— pe/A,  (  *?{t-z/c)  +  Mf{t+s/c)) 

The  density  of  air  is  p,  e  is  the  speed  of  sound  in  air  and  A,  is  the  cross-sectional  area  of  the 
acoustic  tube,  see  Fig.  6a.  Assuming  that  all  the  tubes  are  of  equal  length,  L,  at  the  boundary 
between  tubes  i  and  i+  1,  a  continuous  wave  propagation  is  required. 

,(0,<) 

Using  the  boundary  conditions,  the  transmitted  wave  u£i  and  the  reflected  wave  uf  across  the 
boundary  are  determined. 


«&i(0  -(!+*«)  *+(M  +  *  »f+i(0  (14) 

«i1<+ »)  —  -  *  +  ( i  -*<)•»<+ 1  (0 

Here  f«  L/e  is  the  propagation  time  through  the  tnbe  section  and  kf  is  the  wave  reflection 
coefficient  at  the  junction  of  A(  and  Ai+ ,. 

1 -At)  /  (Ay+ ,+  Af)  (15) 

Since  the  cross  sectional  areas  are  aO  positive,  -l<i<<l.  The  wave  propagation  due  to  the 

discontinuity  in  cross  sectional  area,  is  shown  in  Fig.  6a,b. 

The  lattice  liter  structure  is  obtained  by  normalising  variables  and  grouping  time  delays. 
By  modifying  (14),  the  waves  in  the  i-th  physical  section  at  the  boundary  with  the  (»+ 1  )-tk  sec¬ 
tion  can  be  written  in  terms  of  the  (i+  l)-i*  section. 

u,+(l-r)  -  (  u£,(l)  -  k,  niVt(l) )  /  ( 1  +  h, )  (16) 

urtf+»)-(u*1(l)-*,«4i(0)/(l  +  kt) 

Am  absolute  time  reference  is  sstabBshsd  at  the  output  of  the  last  tube  section,  which  physically 


would  be  at  the  lips.  Assuming  that  the  vocal  tract  model  has  p  tube  sections,  the  time  delay 
from  the  beginning  of  the  i-tk  tube  to  the  lips  is  f,  *  (  p+  1  -»)  r,  so  the  time  variable  in  the 
equation  for  the  i-th  section  is  replaced  by  <  -  f<.  A  scale  factor  c,  is  introduced  to  combine  the 
(  I  +  kj)  factors  from  the  i-tk  tube  to  the  last  (p-th)  section  (lips). 

e>  “  1  +  *,)"■(!+ 

>=»• 

The  lattice  equations  are  obtained  from  (Id)  by  using  the  absolute  time  reference  and  defining 
new  variables,  see  Fig.  6c. 

/,(<)  -  e,  «i+(<  -'-<<) 

MO  —  *iX*  +  *-<<) 

/.(0-A»(0-*i  Mi(<  -2r)  (17) 

MO -MiO -*)-*.  A+»(0 

This  is  the  same  equation  as  that  developed  earlier,  (4)  and  (7),  from  orthogonality  conditions 
except  the  unit  delay  is  2 r  and  the  lattice  sections  are  numbered  in  decreasing  order. 

Although  modeling  of  the  entire  vocal  tract  includes  other  ialuences  due  to  the  vocal  chords 
(glottis)  and  lip  radiation,  the  wave  propagation  in  the  mouth  ideally  follows  the  lattice  structure 
equations.  Studies  have  indicated  that  every  reasonable  vocal  tract  shape  could  be  generated  by 
a  lattice  filter  and  that  the  reflection  coefficients  are  directly  related  to  the  cross  sectional  area  of 
the  vocal  tract  [Markei  and  Gray,  1876], 

For  other  types  of  signals,  if  they  are  generated  by  or  can  be  modeled  as  wave  propagation 
in  a  stratified  medium,  the  lattice  structure  is  intuitively  motivated.  For  physically  generated 
processes,  the  process  is  often  nonstationary  but  there  is  a  limit  to  the  rate  at  which  a  process  can 
change.  The  shape  of  the  vocal  tract  (excluding  the  tips)  can  only  change  at  a  moderately  slow 
rate  determined  primarily  by  the  muscles  in  the  tongue.  Except  when  a  sudden  opening  of  the 
lips  occurs,  the  cross  sectional  area  of  the  vocal  tract  changes  slowly  and  hence  the  reflection 
coefficients  also  change  slowly.  This  slow  time  evolution  can  be  used  advantageously  in  adaptive 
estimation  or  parameter  quantisation. 

An  intuitive  understanding  of  the  significance  of  reflection  coefficient  values  is  possible 


because  the  lattice  filter  structure  can  be  thought  of  as  wave  propagation  in  the  acoustical  tube. 
The  equivalent  tapped  delay  line  coefficients  are  not  as  easily  interpreted.  While  the  reflection 
coefficients  are  limited  to  -l<i,<l,  the  equivalent  tapped  delay  line  coefficients  are  often  a  fac¬ 
tor  of  ten  times  larger.  When  the  reflection  coefficient  b  tero,  the  signal  propagates  without 
change  since  adjoining  sections  in  the  tube  would  have  the  same  cross  sectional  area.  When  the 
reflection  coefficient  is  near  -1,  the  signal  is  apt  to  have  highly  resonant  or  oscillator  characteris¬ 
tics  since  if  the  next  tube  section  is  completely  closed,  ie.  tero  cross  section,  then  the  wave  is 
totally  reflected.  Conversely,  when  the  reflection  coefficient  is  near  -I- 1,  a  decaying  signal  ampli¬ 
tude  is  usually  found  since  if  the  cross  sectional  area  increases  greatly  across  a  boundary,  then 
there  is  full  forward  radiation.  This  connection  between  the  physical  properties  of  wave  propaga¬ 
tion  in  the  acoustical  tube  and  the  analogous  lattice  filter  structure  greatly  aids  in  an  intuitive 
understanding  of  the  effect  of  reflection  coefficient  values  on  signal  characteristics. 


Cross-sectional  Area 


* 
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4.  SAMPLE  DATA  ESTIMATES  OF  REFLECTION  COEFFICIENTS 


The  real  advantage  of  the  lattice  filter  is  for  adaptive  filtering,  where  the  characteristics  of 
an  unknown  process  are  to  be  determined  from  the  observed  data  samples.  The  remainder  of  this 
overview  presents  the  adaptive  lattice  filter  where  each  new  data  sample  is  used  to  update  the 
reflection  coefficients.  In  this  section,  approximation  techniques  that  estimate  the  reflection 
coefficients  based  on  gradient  approaches  or  sample  (block)  data  estimates  of  statistical  quantities 
are  presented.  The  development  of  a  recursive  exact  solution  for  least  squares  estimation,  which 
naturally  produces  a  lattice  structure  filter,  begins  in  Section  5.  A  simpler  set  of  recursive  equa¬ 
tions  using  normalised  variables  is  presented  in  Section  7. 

The  earliest  techniques  for  estimating  reflection  coefficients  assumed  that  the  signal  was 
locally  stationary.  Therefore  sample  data  approximations  were  used  for  the  statistical  definition 
of  the  reflection  coefficients,  (6)  and  (8).  When  the  process  *(■)  is  stationary  with  known  auto¬ 
correlation  function,  the  forward  and  backward  prediction  error  energies  at  each  stage  are  identi¬ 
cal,  (<r/  a j).  Thus  the  reflection  coefficient  for  the  forward  and  backward  predictor  are  the 
same  and  the  lattice  Alter  stage  requires  a  single  parameter. 

MO  -//(<)-*/♦  1  MM  (18) 

*y»(0-M  i /i(0 

The  block  data  techniques  use  a  time  sequence  of  data  and  determine  a  single  prediction 
Alter  for  this  entire  Mock  of  data.  A  single  reflection  coefficient  per  lattice  stage  is  calculated  by 
combining  sample  data  estimates  of  k1  and  k*.  If  the  geometric  mean  of  k1  and  k*  is  used,  then 
the  reflection  coefficient  becomes  the  correlation  coefficient  between  fj  and  t,.  This  parameter, 
k1  wan  originally  called  a  partial  correlation  (PAR COR)  coefficient  [Itakura  and  Saito,  1068).  It 
in  the  normalised  conditional  correlation  coefficient  between  s(<)  and  x(f-;-l)  given  the  interven¬ 
ing  data  mm  pies,  *(l-l) s(l-j'). 


jm-kiUrt 


•-  \ 


The  expression  for  M-i  that  minimizes  E(//+i(0)  +  £(k/+i(0)  “  the  harmonic  mean  of 
k1  and  i*.  This  estimate,  kB  is  computationally  simpler  and  is  related  to  Burg’s  maximum 
entropy  method  (Burg,  1975]. 


£/,(*)  M‘-i) 

tf+i  -  - ^ -  (20) 

1/2  E(/A0+  */M)) 

i=i 

These  two  definitions  for  reflection  coefficients  use  blocks  of  T  data  samples  and  thus  require 


many  computations.  For  kJ+ t,  the  T  samples  of  f  t  and  tj  are  required.  To  calculate  ti+  2,  the 
lattice  filtering  at  stage  j+ 1  must  be  performed  to  obtain  T  new  values  for  /y+ ,  and  Mi- 
Clearly  this  requires  P  filtering  steps  of  2T  samples  to  determine  P  reflection  coefficients. 

Adaptive  gradient  algorithms  for  determining  the  reflection  coefficients  greatly  reduce  the 
computational  complexity  of  the  estimation  technique.  Only  the  prediction  errors  at  the  previous 
time  instant  is  needed  for  the  gradient  methods.  The  block  data  approach  required  all  the  past 
error  values.  Several  techniques  have  been  proposed  for  adapting  the  reflection  coefficients  for 
every  newly  observed  data  sample.  These  techniques  do  not  minimize  any  criterion  but  try  to 
change  the  reflection  coefficient  in  the  direction  of  decreasing  prediction  error  energy  (gradient 
descent).  Two  classes  of  gradient  techniques  either  approximate  the  reflection  definition  of  (20)  as 
the  current  reflection  coefficient  plus  a  correction  or  approximate  the  numerator  and  denominator 
separately.  The  simplest  update  of  the  reflection  coefficient  [Griffiths,  1977]  uses  the  forward  and 


backward  prediction  errors  weighted  by  the  constant  a. 

kj(t+ 1)  -  MO  +  « ( MO  Mi(*-0  +  //-i(0  MO }  (21) 

This  estimate  can  be  improved  by  replacing  the  weighting  factor  by  an  energy  normalized  term, 
1  /  MO  where  MO  1*  accumulated  average  of  and  t£.j(f-l)  [Griffiths,  1978, 

Makhoul,  1978(2)]. 

M*+ 1)  -  MO  +  ( 1  /  *i )  { MO  Mi(‘-0  +  fj- 1(0  MO }  (22) 

MO  -  (i-0  M'-O  +  fi- »(0  +  *A0-0 

Another  adaptive  estimate  [Makhoul  and  Viswanathan,  1978(1),  Makhoul,  1978(2)]  approxi¬ 
mates  the  numerator  and  denominator  of  (20)  separately.  The  same  weighting  factor  a  is  used  for 
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both  terms.  The  ratio  of  these  two  terms  becomes  the  estimate  of  the  reflection  coefficient. 

Cj(t  +  1)  -  (  1  -  a  )«,(<)+  2  /y_i(<)  6y_i(<-l) 

4('+ 1)  -  ( i  -  a )  4(0  +  tM)  + 

<23> 

This  ratio  is  a  biased  estimator  since  in  general  E(x  /  y)y^E(x)/ E(y)  but  simulations  indicate  that 
this  bias  is  generally  very  small  [Honig  and  Messerschmitt,  1981]. 

The  convergence  of  the  lattice  filter  is  much  faster  than  that  of  the  adaptive  tapped  delay 
line  filter  (Satorius  and  Alexander,  1979(2),  Horvath].  This  is  because  the  lattice  filter  tries  to 
orthogonalize  the  input  signal  so  that  the  coefficient  estimates  are  decoupled.  In  fact,  the  conver¬ 
gence  time  is  almost  independent  of  the  eigenvalue  spread  of  the  signal,  i.e.  independent  of  the 
signal’s  spectral  dynamic  range  (Griffiths,  1977],  Quantitative  characterizations  of  the  conver¬ 
gence  properties  of  the  gradient  reflection  coefficient  estimators  (22)  and  (23)  have  been  studied 
(Honig  and  Messerschmitt,  1981].  A  two  stage  gradient  lattice  algorithm  was  compared  with  a 
two  stage  LMS  gradient  tapped  delay  line  filter  to  show  that  it  is  possible  but  unlikely  for  the 
tapped  delay  line  filter  to  converge  faster  than  the  lattice  filter.  A  comparison  of  lattice  estima¬ 
tion  techniques  using  the  gradient  and  block  data  reflection  coefficient  definitions  (21),  (22),  and 
(23)  has  been  presented  in  (Gibson  and  Haykin]. 

This  orthogonalizing  and  decoupling  property  of  the  lattice  is  only  asymptotically  obtained 
using  the  gradient  estimation  techniques  of  this  section.  The  recursive  least  squares  lattice, 
developed  in  the  next  section,  exactly  solves  the  orthogonalization  for  every  new  data  sample. 
The  optimal  solution  is  similar  to  the  energy  normalized  gradient  lattice  (22)  except  that  the 
optimum  weighting  factor  is  computed  (instead  of  the  constant  ff).  Thb  least  squares  lattice  esti¬ 
mation  technique  has  even  faster  convergence  than  the  gradient  lattice  methods  of  above.  How¬ 
ever  as  the  number  of  data  samples  (from  a  stationary  process)  gets  large,  the  results  from  the 
gradient  lattice  and  the  least  squares  lattice  become  similar. 


5.  RECURSIVE  LEAST  SQUARES  LATTICE  ALGORITHM 


The  recursive  Least  Squares  Lattice  algorithm  (LSL)  allows  the  exact  solution  to  the  least 
squares  problem  to  be  updated  for  every  newly  observed  data  sample.  This  adaptive  estimation 
technique  uses  the  properties  of  the  lattice  filter  to  efficiently  implement  the  adaptation.  The 
LSL  algorithm  looks  similar  to  the  gradient  techniques  of  the  previous  section  except  that  optimal 
weighting  factor  are  calculated.  The  LSL  algorithm  is  developed  in  this  section  in  the  context  of 
an  extension  to  the  Levinson  algorithm  for  solving  the  normal  equation. 

The  least  squares  solution  to  a  linear  modeling  problem  can  be  reduced  to  a  simple  set  of 
linear  equations  called  the  normal  or  Yule-Walker  equations.  These  equations,  which  involve  the 
inversion  of  a  covariance  matrix,  has  been  widely  studied  to  reduce  the  computational  burden, 
guarantee  stable  models,  and  handle  nonstationary  processes.  The  linear  predictor  form  of  the 
linear  modeling  problem  is  presented  here. 

The  linear  prediction  model  assumes  that  a  data  sample  at  time  T,  xT  can  be  approximated 
as  if,  a  weighted  sum  of  previous  data  samples.  For  an  p-th  order  linear  predictor  with 
coefficients  (at,  .  .  . ,  a,)  : 

Xf,T  "  "  *1  Xf-l  -  *  —  Op  ZT-p  (24) 

The  coefficients  are  to  be  chosen  so  as  to  minimize  the  mean  squared  error  between  Jj-  and  the 

estimate,  xt,r-  The  p-th  order  covariance  matrix  of  the  process  x()  is  Rp  and  is  composed  of 

elements  r(j. 

Rp  -■  E|X|T:T-,|  *|T:T-p|  1  (25) 

*|T:T-p|  *=  I*T>  *T-1 . zT-p]T 

ru  «■  E[xT.i  xT.j\,  0<  *,j  <  p 

Minimising  the  square  of  the  prediction  error  with  respect  to  the  predictor  coefficients,  {«*} 
requires  that  the  predictor  coefficients  satisfy  (26),  called  the  normal  equation,  where  ot  is  the 


-24- 


minunum  error. 


{(,r  -  *,.r  )’>  -  R,|. 


Unless  the  process  is  a  deterministic  one,  a  unique  solution  to  (26)  exists. 

In  general,  for  a  p-th  order  linear  model,  the  solution  of  the  normal  equation  involves  the 
inversion  of  the  p  by  p  covariance  matrix.  Standard  matrix  inversion  methods  such  as  Gaussian 
elimination  require  0(ps)  computations  (multiplications).  However  for  stationary  random 
processes,  the  covariance  matrix  is  a  Toepliti  matrix. 

rU  -  n  i-j  |.  «<  i,j  <  p 

Using  the  Levinson  algorithm,  the  normal  equation  in  Toeplits  form  can  be  solved  in  0(pz)  com¬ 
putations.  The  Levinson  algorithm  is  an  order  recursive  technique  that  uses  the  solution  for  an 
i-tk  order  predictor  to  generate  the  solution  for  the  (i+  1  )-th  order  predictor.  This  algorithm 
performs  an  orthogonalization  as  discussed  is  Section  3.  The  reflection  coefficients  are  related  to 
the  predictor  coefficients  and  are  generated  as  a  byproduct  of  this  algorithm.  In  the  Toeplits 
case,  the  reflection  coefficients  can  be  determined  directly  without  using  predictor  coefficients 
(LeRoux  and  Gueguen],  This  is  an  application  of  the  Schur  algorithm.  If  the  covariance 
sequence  (r0  ,  rt ,  ,  rf)  is  fed  into  a  growing  order  lattice  filter,  the  reflection  coefficient  at 

each  stage  can  be  determined  by  dividing  the  forward  error  by  backward  error  at  the  input  of 
that  stage  (Morf  et  al,  1977]. 


S.1  Formulation  of  Recursive  Estimates 

In  the  development  of  recursive  least  squares  lattice  algorithms,  two  aspects  of  the  solution 
of  the  normal  equation  are  important.  The  first  aspect  is  the  efficient  inversion  of  the  covariance 
matrix,  Toeplits  and  non-Toeplitz,  that  gives  rise  to  the  order-update  recursions.  Secondly,  the 
time-update  structure  allows  exact  least  squares  solutions  to  be  computed  in  a  recursive  manner 
for  each  new  data  sample.  Thin  enables  the  lattice  algorithms  to  achieve  extremely  fast 


sv-ra-v*:. 


,  •  ,  '  «/  s,' 


v-  .**  . 
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convergence  and  excellent  tracking  capability.  The  first  derivation  of  the  least  squares  lattice 
came  from  extending  the  Levinson  approach  and  was  based  on  the  predictor  coefficients.  Once 
the  recursions  were  obtained,  it  was  noted  that  they  could  be  written  compactly  using  just  the 
lattice  filter  parameters.  A  subsequent  direct  derivation  of  the  recursive  equation,  using  a 
geometric  approach  solely  in  terms  of  reflection  coefficients  and  lattice  parameters  was  reported  in 
[Lee,  1080,  Lee  et  al,  1081,  Shensa,  1081).  An  overview  of  the  first  derivation  is  presented  here 
since  it  gives  insight  into  the  nature  of  the  recursions.  The  development  of  the  order-  and  time- 
update  recursions  for  the  lattice  algorithms  based  on  the  special  structures  of  the  normal  equation 
follows  the  approach  in  [Lee,  1080|. 

First  the  structural  properties  of  the  sample  covariance  matrix  must  be  exploited  as  the 
order  of  the  predictor  changes  and  as  new  time  samples  are  added.  The  covariance  matrix  of 
order  p  for  data  samples  {*y,  j “0, T }  used  here  is  the  prewindowed  form. 


R,.T  -  X,.T  X,.T 


The  covariance  matrix  R|,t  has  the  following  structural  properties.  As  a  new  data  sample 
Zj>i  is  included,  the  new  covariance  matrix  is  composed  of  the  previous  covariance  matrix  plus  a 
matrix  of  special  form.  The  time-update  matrix  equation  is  given  in  (20). 


R<.r  + 


[*r+  v  •  •  • »  1 1 


*T-i+ 1 


Also,  the  covariance  matrix  of  order  i+  1  contains  the  covariance  matrix  of  order  t.  The 
order-update  matrix  equation  for  the  covariance  are  given  in  (30)  where  *  denotes  unspecified  ele¬ 
ments  along  the  outer  row  and  column. 


»•*«  -  [: -  [v :  1 


The  development  of  a  recursive  solution  requires  that  the  normal  equation,  (26)  be  extended 
with  auxiliary  vectors.  The  forward  and  backward  predictors  vectors  are  A,j  and  Cjj  respec¬ 
tively.  A  vector  Dir*9  defined  to  account  for  new  data  samples.  The  matrix  equation  (31),  an 
extension  of  (26)  defines  these  vectors. 

vf.T  I  0  I  *r 

0  |  0  |  xT.x 

Ri.r  I-A. «, r»  Ci,T<  D,  rl  “  I  I  (31) 

0  l  0  |  zr.1+1 

0  I  ffi.T  I  xT-i 

A.-.r  =  |  1.  *i.r  > •  •  •  /  *i,r  I 

Cj,r  *■  I  et,T  ,  •  ■  ■  >  ci ,r>  1 1T 

The  forward  prediction  error,  /,  r  and  the  backward  prediction  error,  bi  T,  are  defined  as  in  (2) 
and  (3)  with  X|j-  j>_,|  defined  as  in  (25). 

Ji.T  *  A<.r  *|rr-.|  (32) 

,  _  T 

°i.T  “  Ci,r*|T:r-il 

To  account  for  the  end  of  the  sample  data  set,  an  auxiliary  vector  D  ,  T  and  related  scalar  7,r 
are  introduced. 

Di.r  “  Ri.T  x\T:T->\  (33) 

r  r  -i  .  . 

7.-,r  —  D*,r  *|rr-»i  —  *|rr-<|  «.\r  z|rr-.|  (34) 

This  parameter  7ir  can  be  interpreted  as  a  likelihood  variable  and  is  limited  to  the  range 
0  <  7i  r  <  1  ,  see  Section  5.5. 

5.3  Order  Update  Equations 

As  in  the  Levinson  solution,  an  efficient  means  of  determining  the  n-th  order  solution  is  to 
develop  recursive  equations  for  updating  the  predictor  order  (at  a  fixed  time).  Following  the 
usual  development  of  recursive  equations,  assuming  the  predictor  vectors  A ,  r  and  C,  r  are 
known,  the  predictors  of  order  i  +  1  are  to  be  determined.  Thus  the  vectors  A,+  1>r  and  Ci+ 
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must  satisfy  the  following  normal  equation. 


®i+i,r  [Aj+iji 


•+i.r  I  0 

0  I  0 

I 

6  |  6 

0  I  *i+i,r 


Since  A ,+  jt j  b  to  be  assembled  from  A(y  and  C i  T.u  these  i-tk  order  predictors  are  aug¬ 
mented  with  a  final  (or  initial)  zero  to  make  an  »+  1  vector.  From  (30),  the  augmented  A ,  r 
satbfies  the  new  normal  equation  except  for  the  last  entry. 


•>.r  {*"o  r]  -  fV  :]  b 


I  Oi,T>  0  >  •  •  •  /  0.  A,+  lir] 


A(+i,r  ”  (f«*<  row  of  R , 


•iH  [V] 


Similarly,  augmenting  C^r-i  with  a  leading  zero  satisfies  the  normal  equations  except  for 


the  first  entry. 


R<+l-r  W  “  ['  R  *r-»]  [c°r.,]  -  lr'+l'r'  0 .  °>  ’ 

r,+  i.r  -  l [/*>*»  row  of  Rl+i.rl  |c,°r_i] 


Since  ,  r  is  symmetric,  it  can  be  seen  that  Af+i,r  “  r/-n,r  IBu75). 

By  appropriately  combining  these  two  equations  to  cancel  the  leading  or  final  term  in  the 
normal  equation,  the  order-update  equations  are  obtained.  Multiplying  (37)  by  A<+i.r  / 
and  then  subtracting  the  result  from  (35),  the  order-update  recursion  for  A,  r  and  <riiT  are 


obtained. 


*-r  “  [V]  ~  ****  I 


ffLl,T  "  *!,T  -  &hl,r  /  ai,T-l 


.-.V 
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Similarly  the  order-update  recursion  for  C ,  j  and  <r,  T  are  obtained. 


Ci+i.r  -  Jc ®r_,j  "  [V]-,  /  ^ 

(40) 

»  *  A2  ,  l 

“  °i.T- 1  -  &i+l,T  /  ai,T 

(41) 

When  the  order-update  equations  of  the  predictors  are  premultiplied  by  [xy,  .. 

the  1st- 

tice  equations  result.  The  reflection  coefficients  are  i/+  ,  r  and  £,>  17. 

fi+ 1. r  —  fi,r  -  *i+i.r4i.r-» 

(42) 

*i+l ,T  “  ”  *<+  l,r  fi.T 

(43) 

J  A  ,  / 

*1+  i.r  “ 

(44) 

*i+i.r  ™  A»+  i  t  1  ^j,r-i 

(45) 

So  far  the  development  haa  been  parallel  to  the  non-adaptive  approach  in  Section  3.  To  develop 
adaptive  solutions,  time-update  equations  for  the  predictors  must  be  developed.  Before  proceed¬ 
ing  to  time-updates,  the  order-update  recursion  for  D ,  r  is  developed  in  a  similar  fashion. 


The  last  element  of  D  ,>  lt  7  can  be  found  from  the  last  row  of  R 1, 
lost  element  D  i  t  ■■  lost  element  (  R  »>  1,  r  )  *  |  r.  r-«-i| 


*r-o  *r 

T,  by  (32)  and  (33). 
“  &«+ 1  ,r  /  1+ 1, r 


Since  the  last  element  of  C  l+  tr  is  1  and  the  last  element  of  D  (+  }  T  has  been  determined,  the 
order-update  equation  for  D  i+  i  r  becomes  (46). 


D<+i,r  “  f*r]  +  J(+i  r/rf+Iir 

The  order-update  for  li  t  is  determined  by  premultiplying  (46)  by  [xT . 


7<+i,r 


li.T  +  *<+»,r  /  <7<+i,r 


(46) 


(47) 


IJ  Time  Update  Equations 

Next  the  time-update  of  the  covariaace  matrix  is  used  to  determine  time-updates  for  the 
predictor  vectors.  From  the  time-update  of  the  covariance  matrix  (29)  and  the  definition  of  the 


1 

M 


forward  prediction  error  (32); 


/  r  _ 
viT  *r+ 1 

0  *r  r 

Ri,r+i^i,r  **  +  [*t  >  ■  •  •  >  *r-«|  A,r 

0  *r-« 

, 0  pr-i+i 

The  auxiliary  vector  D,  ,r  is  used  to  account  for  the  new  data  samples. 


[o.°i.rJ  “  [*  R,*.r]|b  °lir] 


rr< 

!*r-.+ 1 


The  time-update  recursion  for  A  ,r  is  obtained  from  (48)  and  (49). 


A,,r+i  “  A 


i.T  ~  |p,  t  f|  [XT  >  •  •  •  »  *T-i!  Aj  ] 


By  premultiplying  (50)  by  [zr+ 1.  •••,  *r-«>  lj  and  using  the  definition  of  7,.i  r,  (34)  the  expres¬ 
sion  can  be  simplified. 

\xT  , . . . ,  *r-«|rA<,r  ■«  /.,r+i  /  (1  —  7.-i,r)  (51) 

The  time-update  for  Af  j  becomes  a  simple  expression. 

A,.fi  -  A,.t  -  [D|°ir]  7T^7  <“) 

From  the  preceding  relation,  the  time-update  for  ajj  can  be  determined  using  (52)  and  (29). 

ff!,T+i  "  Ajir+lJ?jir+idjj+, 

/  /  /<?r+ 1 

fi.r+i  *i,r  +  - - -  (53) 

1  -  7<-i,r 

By  applying  similar  techniques,  the  time-update  recursions  for  C  i  T  and  a*  T  are  determined. 


,r..  "  C‘’  -  [  0  J  1  - 


»  » 


ffi,T+l  ™  *i.T  + 


1  -  7<-i,r+i 


To  recnrsively  update  the  reflection  coefficients,  the  time-update  equation  for  A  is  needed. 


10  C,>r)R,+  1<r+1|4,br+1J  A<+ir+j 


I"  » i 


.•  r. 
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By  using  the  covariance  time-update  (29),  the  time-update  for  C  ,r  and  the  forward  and  back 
ward  predictors,  the  time-update  equation  is  obtained. 

*i.r  fi.T+  1 


A<+i,r+i  ™  Aj+t,r  + 


(57) 


1  -  7m.  r 

The  exact  time-update  of  ^  r  is  a  time  average  of  the  cross-correlations  between  6,  r  and 


fi'T+v  with  the  special  gain  factor 


1 


This  relates  Al+  XT  to  the  partial  correlations 


1  -  7i-i.  r 

discussed  in  the  previous  section. 

The  development  of  the  time-update  equations  allows  the  influence  of  all  prior  data  to  be 
accumulated  into  the  current  parameter  estimates.  If  the  process  that  is  being  approximated  con¬ 
tains  time  varying  attributes,  then  it  is  desirable  to  weight  more  recent  observations  more 
heavily.  An  exponential  weighting  factor  X  on  the  accumulated  covariances  can  be  included  in 
the  development  of  this  algorithm.  Typical  values  of  X  are  from  .98  to  1.00  (corresponding  to  full 
weighting  of  past  samples).  The  algorithms  presented  in  subsequent  sections  include  this 
exponential  weighting  factor. 

8.4  Exact  Least  Squares  Lattice  Recursions 

The  recursions  so  far  have  focused  on  the  predictor  vectors,  Ajfj,  and  D|,t-  For  a  P- 
th  order  prediction  filter,  these  recursions  require  0(P*)  operations  per  time  sample  since  all  the 
predictor  coefficients  change  in  an  order-update.  However  for  the  lattice  structure,  due  to  its 
orthogonalising  nature,  only  the  i-th  reflection  coefficient  changes  in  the  i-th  order-update.  The 
exact  least  squares  recursion  can  be  written  directly  in  terms  of  lattice  Biter  variables  which 
require  only  O(P)  operations  to  update  per  time  sample. 


The  order-apdate  recursions  for  the  lattice  filter  variables,  f(j  and  bi  r  was  developed  in 
(42-45).  The  order-updates  for  and  <t*t  are  given  in  (39)  and  (41)  and  the  time-updates  are 
given  in  (53)  aad  (55).  The  reflection  coefficients  defined  above  depend  on  A^t.r-t-1  which  is 
related  to  partial  correlations.  Here,  the  time-update  for  Af+i,r+i  ()  »  required  to  augment  the 
correlation  for  the  new  data  sample.  These  updates  also  require  the  order-update  of  7<ir,  deter¬ 
mined  in  (47). 
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These  recursions  are  seen  to  compute  the  sample  cross-covariance  of  the  forward  and  back¬ 
ward  prediction  errors,  using  the  optimal  weighting  1/(1  -  Tf,_i  r  ).  From  Section  4,  the  gradient 
lattice  equations  have  the  same  form  as  above  except  that  the  exact  recursive  least  squares  solu¬ 
tion  calculates  and  uses  the  optimal  weighting  factor  for  the  new  data  sample. 

The  complete  set  of  order-update  and  time-update  recursions  to  obtain  the  exact  least 
squares  lattice  predictor  (LSL)  is  presented  in  ALGORITHM  3.  When  starting  the  lattice  filter, 
only  the  stages  that  receive  data  are  executed  until  P  data  samples  have  been  observed,  ie. 
min(T,P)  filter  stages  are  used  where  T  is  the  data  sample  number. 


1  - ^Vy*i**^ > » *•'. '»  .*•  »'*  ^  • 


■vv.“ 
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Input  Parameters 


Variables 


Initialisation 


ALGORITHMS: 

Recursive  Least  Squares  Lattice  (LSL) 
Scalar  case  for  exponentially  weighted  data 


P  »  maximum  order  of  lattice 
X  a*  exponential  weighting  factor  (usually  .98  to  1.0) 
xT  ss  data  sample  at  time  T 

A|,r  ”  partial  autocorrelation  coefficients 
7 i'T  *  likelihood  variable 

/i,r  .  Vr  “  forward  (backward)  prediction  errors 

t  >  ffi,T  =  forward  (backward)  prediction  error  covariances 
kjT  ,  ki  t  “  forward  (backward)  reflection  coefficient 

/i  /  »  2 

JOfi  “  #8.0  “  *0  t  “  <*0.0  “  *0 

A,.,  -  0,  7-m  =0,  1<  i  <P 


ITERATION  FOR  EVERY  NEW  DATA  SAMPLE 

For  data  sample  xr  and  previous  results  A,  7,  4  ,  a  ,  <r* 

/o.r  "  *o.r  “  *r 
/  •  *  /  2 
^o.r  ™  tfo.r  ™  a  <ro,r-i  +  *r 

For  each  stage  of  the  lattice,  i  —  0  ,  mm(r,P)-l 

A,+  i,r  “  *  Aj+1,r-i  +  *»,r-i /«.r  /  (1  -  7«-i.r-i) 
7<,r  "  +  *«?r  /  *<,r 


*/+i.r 

- 

^i+i.r 

.  t 

1  ffi.T 

- 

A<+  i,r 

,  * 

/<+i.r 

- 

fi.T  ~ 

*#+i,r 

- 

h,r-i  ■ 

■  */+i ,r/<(r 

wAen  T  <  P  aj+  i  r  —  <r*r  -  kt+i,r  A<+i,r 

*  *  J  A 

#k#  ff!+  i,r  ■■  ^  +  /<+i,r/(l  -  7<,r-i ) 

*<+i,r  ™  *  <7<+i,r-i  +  *<+  i,  /  /  ( 1  -  7<,r  ) 


wv 


ma 


S.S  A  Likelihood  Variable 

Information  about  changes  in  the  nature  of  the  observed  process  can  be  determined  from 
lr,T-  This  variable  can  be  interpreted  as  a  sample  data  approximation  of  a  statistical  likelihood 
variable.  For  a  zero  mean  Gaussian  random  process,  z,  the  joint  distribution  for 
{iT,  xr_i,  ....  if-p)  is  given  by  (58)  with  the  covariance  matrix  as  defined  in  (25). 

1  T  ”1 

p(*T>  *r-r)  “  |2xRp|*l/,J  exp  {  -  —  X|x.t-P|  Rp  *|T:T-p|  }  (58) 

The  determinant  of  the  covariance  matrix,  |RP  |  is  related  to  the  reflection  coefficients  and  the 
process  variance  R0  (Market  and  Gray,  1976]. 

|Rp|  “flofifl-  tf)  (59) 

1  =  1 

The  logarithm  of  (58)  becomes  a  log-likelihood  function  composed  of  two  parts.  The  first  two 
terms  depend  on  the  covariance  of  the  process  and  the  third  term  relies  on  the  observed  data  sam- 


P  T  -1 

log— likelihood  *■  {  In  R o  +  J]  In  (  1  -  K*  )  }  +  *|T:T-P|  Bp  *|T:T-P|  (69) 


The  variable  7P  T  obtained  in  the  exact  least-squares  recursions  can  be  interpreted  as  the  sample 
estimate  of  the  third  term  in  (60).  The  definition  of  7p  p  uses  the  sample  estimate  of  the  covari¬ 
ance  matrix,  Ap,r  instead  of  the  known  covariance  matrix,  Rp.  Thus  7p  r  is  a  measure  of  the 
likelihood  that  the  P  most  recent  data  samples,  (zr,  .  .  .  ,zr_p)  come  from  a  Gaussian  process 
with  sample  covariance  Rp  T  determined  from  all  of  the  past  observations  (z,,  0 <j<T}.  Since 
0  <  7 p  T  <  1  ,  a  small  value  of  7 p  T  indicates  that  the  recent  data  samples  are  likely  observa¬ 
tions  from  a  Gaussian  process  with  covariance  Rp_f.  However,  a  value  of  7 p  T  near  one  implied 
that  given  the  current  Gaussian  process  assumption,  the  observations  are  unexpected;  either  the 
new  observations  come  from  a  different  Gaussian  process  due  to  a  time  varying  nature  of  the  phy¬ 
sical  process  or  there  is  a  non-Gaussian  component  in  the  observations.  So  7 >  f  can  be  used  as  a 
detection  statistic  for  changes  in  the  process  characteristics  or  for  unexpected  (non-Gaussian) 
components  in  the  observations.  Simulations  indeed  demonstrated  that  7 p  T  does  take  values 


*0 


close  to  one  when  sudden  changes  in  the  observations  occurred.  In  the  LSL  algorithm,  'ip  r  acts 
as  an  optimal  gain  control  since  the  new  observation  influences  the  accumulated  estimate  by  a 
factor  of  (1  -  With  this  gain  factor,  changes  in  the  process  statistics  can  instantane¬ 

ously  influence  the  estimates  more  than  just  being  averaged  with  past  observations.  Simulation 
results  that  demonstrate  this  behavior  on  synthetic  signals  and  speech  signals  are  shown  in  Sec- 


6.  JOINT  PROCESS  LATTICE  FILTER 


Many  practical  problems  require  the  joint  interaction  of  two  or  more  processes  rather  than 
the  prediction  of  a  process  based  on  its  own  past  observations.  For  example,  a  channel  equalizer 
in  its  adaptation  phase  uses  the  received  distorted  signal  from  the  channel  and  the  actual 
transmitted  channel  symbols  to  determine  the  channel  distortion  characteristics.  In  a  noise  can* 
celler,  the  information  signal  plus  noise  and  a  reference  or  noise  signal  is  used  to  extract  the  infor¬ 
mation.  The  estimation  of  autoregressive  moving  average  (ARMA)  processes  with  known  input 
involves  the  joint  interaction  of  two  processes.  iQpgenerak,  a  multi-channel  problem  can  be  formu¬ 
lated  as  a  single  vector  process  problem  which  can  be  solved  by  extending  the  previous  least 
squares  solution  to  the  vector  case.  The  one  channel  prediction  lattice  of  Section  5  is  extended 
here  to  include  a  second  related  channel.  This  joint  process  recursive  least  squares  provides  very 
fast  tracking  or  adaptation  for  channel  equalization  or  noise  cancelling. 

When  one  process  y  is  to  be  estimated  from  observations  of  a  related  process  x,  it  is  possi¬ 
ble  to  combine  them  into  a  joint  process  ( x,y ),  that  can  be  solved  as  a  joint  process  lattice  filter. 
The  exact  least  squares  solution  for  joint  process  estimation  is  an  extension  of  the  development  in 
Section  5.  A  new  prediction  error  is  defined  that  includes  samples  from  both  processes.  The  joint 
prediction  error,  jPT  is  the  error  in  estimating  Vr  from  {xT,  zr_,,  ...,  xT_P}  where  {gf}  are  the 
prediction  coefficients  obtained  by  minimizing  the  sum  of  the  squared  errors. 


ir.T  “  VT  +  S  *T-i  (61) 

M 

The  solution  of  (61)  can  be  formulated  in  terms  of  the  lattice  structure  just  as  the  single  process 
predictor  (32)  was  translated  into  (42).  A  prediction  lattice  filter  (LSL)  for  xj  performs  a  Gram- 
Schmidt  orthogonalization  of  {xr_,}  into  the  mutually  orthogonal  backward  prediction  errors 


The  advantage  of  using  the  orthogonal  { hr-< }  instead  of  {x^}  in  (61)  is  that  the  joint 
predictor  coefficients  {iD  become  decoupled  from  one  another  so  faster  convergence  is  possible. 

The  joint  process  lattice  solution  involves  the  LSL  for  the  x  process  and  a  similar  lattice 
recursion  of  the  joint  prediction  error,  j\t-  From  the  LSL  for  the  x  process,  at  the  i-th  lattice 
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stage,  b,  T  is  the  backward  prediction  error,  a,  j  is  its  variance  and  the  likelihood  variable  is 
'li-i.r-  A  new  cross  correlation  term,  A /r,  similar  to  (57),  can  be  defined  between  signals  avail¬ 
able  after  the  i-tk  lattice  stage,  ji.yr  and  bi  T.  This  new  term  can  be  recursively  updated. 


■ 

A  j  _ 

V  A  *  l  j 

ji-l.T  h.r 

i 

**i,T  = 

1  -  7i-i,r 

The  recursion  for  the  joint  prediction  error  /t,r  is  similar  to  (42)  except  that  quantities  after  the 
i-th  lattice  stage  are  used.  The  initial  condition  is  i  r  =  yT  and  the  output  is  jf,  T. 

Ji.T  “  jif-l.T  -  (  &i,T  /  °i,T  )  K.T  (63) 

The  previous  single  channel  LSL  equations  augmented  with  (62)  and  (63)  form  the  complete  solu¬ 
tion  to  the  joint  estimation  problem  and  lead  to  the  joint  process  lattice  filter,  Fig.  7. 

For  noise  cancelling  problems,  noisy  data  containing  the  signal  of  interest,  (yr)  are 
observed  together  with  the  noise  estimate  or  reference  signal,  {zr},  see  [Satorius  et  al,  1979(1), 
Griffiths,  1979(2)].  For  channel  equalization  problems,  (j/r)  is  a  known  training  sequence  sent 
though  the  channel  and  {zr}  is  the  distorted  channel  output.  Applications  of  the  joint  process 
lattice  estimation  algorithms  to  adaptive  data  equalization  have  been  investigated  (Satorius  and 
Alexander,  1979(2),  Satorius  and  Pack,  1981]. 

The  ARMA  estimation  problem  with  known  input  and  with  bootstrap  estimated  input  was 
formulated  as  a  two  channel  lattice  filter  in  (Lee  et  al,  1982].  For  an  input  process  y,  the  output 
ARMA  process  is  z  generated  in  the  following  manner. 

p  P 

*r  +  2  ai  xr-t  *“  Vt  +  2  Vr-i 

(=i  i=i 

A  prediction  equation  can  be  written  for  the  process  z,  if  the  input  y  is  considered  known.  This 
predictor  follows  from  Section  5  but  now  includes  a  weighted  combination  of  past  inputs  y . 

p  P 

*r  “  -  2  °»  *T-i  +  2  VT-i 

i=i  i=i 

Similarly,  a  prediction  equation  for  the  y  process  with  z  known  is  generated  by  extending  (61)  to 
include  a  weighted  combination  of  past  inputs  y. 

p  P 

$T  ™  -  2C.  VT-i  +2  xT-i 

1=1  i=l 


*  ■  v.£v:v/.v:  - 
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The  vector  process  [xj->Vr]  has  a  structure  similar  to  the  scalar  AR  processes  discussed  earlier. 
Prediction  errors  are  now  vectors  and  covariances  are  matrices.  The  ARMA  lattice  estimation 


algorithm  follows  the  scalar  LSL  or  SQNLSL  algorithm  but  the  quantities  are  now  vectors  and 
matrices.  The  reflection  coefficient  has  become  a  two  by  two  matrix.  Further  details  are  found  in 
[Lee  et  al,  1982). 


Joint  Process  Lattice  Filter 


7.  SQUARE  ROOT  NORMALIZED  LEAST  SQUARES  LATTICE 


The  complexity  of  the  recursions  can  be  reduced  and  the  numerical  properties  of  the  vari¬ 
ables  improved  by  rewriting  the  Least  Squares  Lattice  algorithm  in  terms  of  normalized  variables. 
The  Square  Root  Normalized  Least  Squares  Lattice,  (SQNLSL)  developed  in  [Lee,  1980]  has  only 
three  recursions  per  order  for  each  time  sample  where  all  three  variables  have  unit  variance.  The 
reduction  of  the  LSL  into  the  SQNLSL  requires  two  types  of  normalizations.  A  variance  normali¬ 
zation  scales  the  variables  by  their  respective  variances.  A  normalization  by  the  optimal  weight¬ 
ing  factor  using  the  likelihood  variable,  7  is  also  necessary.  A  brief  development  of  the  SQNLSL 
is  presented  here;  see  [Lee  et  al,  1981]  for  more  details. 

The  forward  and  backward  prediction  errors  when  normalized  become  i/,  T  and  r\i  T  respect¬ 
fully.  The  normalizing  factors  are  the  square  roots  of  the  variances,  <r{  r  and  <r,*r  and  the  square 
root  of  the  optimal  weighting  factor  (  1  -  7,.,  r-1 ). 

vi,T  -  fi.T  (  °l.T  )  '1/2  (  1  -  7.-1,  r-1  )~1/2  (64) 

Vi.r-i  =*  *.,r-i  (  ffi.  t-i  )  ^2  ( l  -  7i-i,  r-i  )-1^2 

By  combining  the  two  reflection  coefficients  from  the  LSL,  the  normalized  partial  correlation,  p,  r 
is  defined  like  a  correlation  coefficient.  This  single  parameter  is  the  new  reflection  coefficient. 

Pi+  i,t  —  (  °!,t  )  ~1^2  A,+  i,r  (  )  1/2  (65) 

First  the  recursion  for  the  partial  correlation  (65)  will  be  developed  from  the  LSL  algorithm. 

The  variances  of  the  prediction  errors  (with  exponential  weighting  X  )  has  a  time-update  recur¬ 
sions  given  by  (53). 

,  /  /  /*r 

*  *  °i,T  ~  - - 

*  -  7i-i,r-i 

By  dhridiag  by  (rjr  and  using  the  definition  for  Vi,r>  the  time-update  for  the  variance  can  be 
related  to  the  new  variables  (66).  A  similar  relation  (67)  exists  for  the  backward  prediction  error 
variance. 

*  I  al,T  ™  1  ~  vf,T  (66) 

-  l-nlr  (67) 

The  time-update  recursion  for  the  normalised  partial  correlations  are  obtained  by  substituting  the 


-40- 


expression  for  Aj+  i  r  (57)  (including  the  exponential  weighting  X  )  in  (65). 


Ai+i.r  ™  Aj+i.r-i  + 


^i,r-l  fi.T 

1  -  7i-i.r 


Pt+i.r  =  (  °!.t  )  *^2  (  ^  Ai+i.r-i  +  7  1  )  I/<2  (68) 

1  -  7i-i.r-i 

The  first  term  is  replaced  by  the  definition  of  Pi+i,r-i  and  second  term  is  vi  r  17,  r_,. 

Pi+i.r  ™  *  (  °!.t  )  ~1/2  (  ) l//z  Pi>i.r-i  (  ai,T-z  )  ’^2  (  ai,T-\ ) 1/2 

+  ui,T  Vi.T-l 

Using  (66)  and  (67),  the  new  time-update  equation  for  the  normalized  partial  correlation  simplifies 


to  (69). 


1.T  —  ( 1  -  I'f.r  )1/2  /><+  i.r-1  ( 1  -  nf,T-i  )1/2  +  ^,r  n.,r-i 


Now  the  lattice  recursions  can  be  written  in  terms  of  these  new  variables.  The  order-update 
recursions  for  the  forward  prediction  errors  (42)  can  be  written  using  the  normalized  partial  corre¬ 
lation. 

(1  -  7i,r-i)1/2  "i+i.r  (*/+i ,r)1/2  *  (*,'r)1/2  Kr  -  Pi+i.r  *7., r-i)  (1  -  7i-i,r-i)1/2  (70) 

To  simplify  this  expression,  two  order-update  equations  from  the  development  of  LSL  are  needed; 

for  the  likelihood  variable  (47)  and  for  the  prediction  error  variances. 

( 1  -  7<,r  )  =  (  1  -  7i-i,r  )  ( 1  -  tf.T  ) 

vf+l.T  /  °!,T  =  1-P?+1,T 

Using  these  relations,  (70)  can  be  reduced  to  a  simple  expression  for  the  normalized  forward  pred¬ 
iction  errors  (71).  A  similar  development  for  the  backwards  prediction  error  leads  to  (72). 

|/<+i,r  *■  ( 1  -  P?+i,r  r1/2  (  vi.T  ~  Pt+i.T  Vt.T-i )  ( 1  —  *7?, r-z  )-1/2  (71) 

Vi+l.T  “  (I  -  Phl.T  )-1/2  (  Vi,T-\  -  Pi+l,T  vi,T  )  (  1  -  vi,T  )"'/2  (72) 

The  lattice  recursions  have  now  become  three  equations,  (69),  (71),  and  (72)  that  compute  the 
normalized  prediction  errors,  {v}  and  {17},  and  the  reflection  coefficients,  {p}  for  each  lattice  stage 
and  for  every  data  sample.  Proper  initialization  is  required  to  start  the  recursions  with  unit  vari¬ 
ance  quantities. 


. --.v, .v.v .v.v.v,’. >,o,.v.v.vv>v  \- v 
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The  reflection  coefficients  in  the  SQNLSL  still  have  magnitudes  bounded  by  one  but  now 
the  prediction  errors  are  also  bounded.  The  complexity  of  the  lattice  recursions  has  been  reduced 
from  six  recursions  to  only  three  recursions  per  order  and  time  update.  Square  root  operations  are 
required.  They  can  be  efficiently  computed  by  bit  recursive  algorithms  such  as  the  CORDIC 
technique  (discussed  in  the  next  section).  Simple  recursions  and  their  potentially  better  numerical 
behavior,  makes  the  SQNLSL  algorithm  preferable  over  the  unnormalized  LSL  version. 

The  SQNLSL  just  developed  applies  for  exponentially  weighted  data.  However  the 
exponential  weight  X  in  (68)  is  not  evident  in  these  recursions.  By  combining  tbe  three  time- 
update  recursions  for  of  T,  af  t  and  A(+l,r  into  one  recursion  for  pi+  lT,  the  effect  of  X  is  carried 
through  unseen.  When  a  new  data  sample  is  used,  the  exponential  weighting  is  applied  to  the 
sample  variance  estimate.  ALGORITHM  4  summarizes  the  Square  Root  Normalized  Least 
Squares  Lattice  (SQNLSL)  estimation  method.  The  sample  variance  RT  is  initialized  to  some 
value  a  to  avoid  dividing  by  zero. 

Although  SQNLSL  is  a  very  powerful  and  compact  algorithm,  the  necessity  of  computing 
square  roots  can  lead  to  problems.  The  fixed  point  error  analysis  of  this  algorithm  [Samson  and 
Reddy]  indicated  that  finite  wordlength  arithmetic  computation  of  the  square  roots  lead  to  small 
biases  in  the  reflection  coefficients.  This  bias  was  more  predominant  than  the  variance  of  the 
error  in  the  estimate  and  generally  quite  small.  The  bias  increased  as  the  wordlength  became 
shorter  or  the  exponential  weighting  factor  X  approached  one. 
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ALGORITHM  4: 


Square  Root  Normalized  Least  Squares  Lattice  (SQNLSL) 
Exponentially  Weighted  Scalar  Data 


Input  Parameters 

P  *■  maximum  order  of  lattice  filter 
X  *  exponential  weighting  factor  (usually  .98  to  1.0) 
a  *  prior  variance 
it  —  data  sample  at  time  T 

Variables 

Rt  99  estimated  variance  of  x 
Pi,  t  —  reflection  coefficients 
ViT  =  normalized  forward  prediction  error 
tfi  T  =®  normalized  backward  prediction  error 

Initialization 

R,  —  a  +  it 
Vt.s  =  V#  =  Ztfy/R, 

Pi. ,  =0  1  <i  <P 

ITERATION  FOR  EVERY  NEW  DATA  SAMPLE 

New  data  sample  xT  and  previous  results  Vi,r-i  >  P>+  i,r-i »  Rt- i 

Rt  —  X  +  x } 

v4.T  *■  t>4,T  “  *t/>/Rt 

For  each  stage  of  the  lattice,  i  ™  0  to  m>n(r,P)-l 

y/1  ~  vi,T  1  Pi+  l.T-l  +  Vi.Tli.T-l 

Vj.T  -  Pi+l,TVi,T-\ 

"  phl.T  n/1  - 


Pi+i.r  “ 
‘'i+i.r  * 


8.  COMPUTATIONAL  COMPLEXITY  AND  CORDIC  ARITHMETIC 


The  complexity  of  the  lattice  filter  is  greater  than  the  equivalent  tapped  delay  line  filter, 
However  the  lattice  filter  has  many  advantageous  properties  not  shared  by  the  simpler  filter. 
Similarly  the  adaptive  lattice  filter  has  a  few  more  operations  than  an  adaptive  tapped  delay  line 
filter.  Table  I  compares  the  computational  complexity  of  several  adaptive  estimation  algorithms. 
The  lattice  methods  require  three  to  six  times  as  many  computations  as  the  simplest  adaptive 
tapped  delay  line  filter  (LMS).  However  this  increase  in  computational  complexity  provides  for 
substantially  faster  convergence,  better  numerical  properties  of  the  coefficients,  and  an  assurance 
of  a  stable  filter.  The  complexity  of  several  adaptive  algorithms  is  presented  in  the  following 
table.  The  scaling  by  a  constant  weighting  factor,  eg.  X  or  0  is  usually  approximated  as  a  shifting 
by  a  power  of  two.  Thus  this  fixed  scaling  is  not  included  in  the  count  of  operations.  The  LMS 
algorithm  is  the  tapped  delay  line  gradient  least  mean  squares  technique.  The  gradient  lattice 
algorithms  is  (22)  and  (18).  ALGORITHM  3  is  denoted  LSL  and  ALGORITHM  4  is  called 
SQNLSL.  The  number  of  operations  are  counted  for  executing  a  single  filter  stage  on  a  single 
data  sample.  To  process  T  data  samples  in  an  N-th  order  filter  would  require  NT  time  as  many 
computations. 


TABLE  1  Computational  Complexity 

Algorithm 

X 

-r 

y/ 

± 

LMS 

2 

0 

0 

2 

Grad.  Latt. 

6 

1 

0 

6 

LSL 

6 

6 

0 

7 

SQNLSL 

10 

2 

3 

6 

The  SQNLSL  algorithm  is  the  most  complex  of  the  algorithms,  requiring  three  square  roots, 
ten  multiplications,  and  two  divisions  to  execute  an  update  for  each  stage  in  the  lattice  for  every 
new  data  sample.  However,  the  SQNLSL  has  a  very  compact  form  with  only  three  equations 
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involving  variables  that  have  constrained  magnitudes.  The  implementation  of  the  SQNLSL  equa¬ 
tions  in  hardware  would  require  special  multiplier  hardware  for  fast  execution  or  would  require 
considerable  execution  time  on  general  purpose  microprocessors.  Even  using  shift  and  add  instead 
of  multiplication  and  Newton’s  method  for  square  root  operations  would  require  considerable  exe¬ 
cution  time. 

By  interpreting  the  lattice  equations  as  rotations,  an  efficient  realization  of  the  square  root 
normalized  lattice  algorithm  using  CORDIC  arithmetic  was  developed  in  [Ahmed,  1981(1)].  The 
Coordinate  Rotation  Digital  Computer  (CORDIC)  developed  in  [Voider]  is  an  iterative  algorithm 
for  computing  trigonometric  functions,  multiplications,  divisions,  and  square  roots.  The  CORDIC 
algorithm  interprets  the  above  functions  as  rotations  of  a  vector  in  different  coordinate  systems. 
The  rotation  is  implemented  by  a  sequence  of  shift  and  add  operations.  This  type  of  arithmetic  is 
not  new;  it  has  been  used  to  compute  trigonometric  functions  and  their  inverses  in  hand  held  cal¬ 
culators.  Many  other  signal  processing  algorithms,  such  as  DFT  and  matrix  inversion,  can  also  be 
implemented  in  arrays  of  CORDIC  processors  (Ahmed  et  al,  1982]. 

8.1  CORDIC  Arithmetic 

The  well  known  equation  for  rotating  a  vector  [z,,  y,]r  to  a  new  vector  [zf+1,  y1+  Jr  uses 
the  tine  and  eorine  of  the  rotation  angle  9. 


Four  multiplications  by  two  trigonometric  quantities  are  required.  This  operation  can  be  made 
more  amenable  to  fast  computer  implementation  by  modifying  this  equation  into  a  sequence  of 
small  rotations  of  a  specific  form,  each  implemented  using  only  additions  and  shifts. 

CORDIC  arithmetic  was  unified  into  a  single  equation  [Waltherj  that  allows  rotations  on 
either  a  circle  (mol),  along  a  line  (n«0),  or  along  a  hyperbola  (mo-1).  The  incremental  unit 
of  rotation  at  the  i~th  iteration  is  the  predetermined  sequence  (6; }  and  (i,  *  ±1  that  determines 
the  direction  of  rotation. 


[*<♦1]  _  f  1  »»M<]  M 
[ft+ij  1  J  [V(J 


(78) 


r\9r9^. 


Pi 
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The  {5,}  must  be  chosen  to  satisfy  certain  constraints  that  assure  convergence  of  the  iterations. 
To  obtain  computational  efficiency  on  current  computer  hardware  using  binary  representations,  a 
negative  integer  power  of  two  is  chosen  for  St.  This  allows  the  multiplication  by  £,  to  be  per¬ 
formed  as  a  shift. 

The  effect  of  (78)  is  a  rotation  and  a  change  of  scale  interpreted  in  the  appropriate  coordi¬ 
nate  space.  The  vector  (x0,  y0]r  can  be  represented  as  a  generalized  radial  component  R0  and  a 
generalized  angular  component,  0O. 

R  o  —  V*o  +  mvT 

*o  «-  \ZrrT  tan ~l(i/oVm  /*o) 

For  the  circular  rotation,  this  is  a  true  polar  coordinate  representation.  Performing  the  operation 
in  (78)  scales  the  radial  component  by  ra-  =  s/l  +  m6f  and  changes  the  angular  component  by 
0,  =»  m~1/2  tan ‘‘(fi.Vm’)  .  After  p  iterations,  the  new  radial  and  angular  components  are  Rr  and 

V 

R,  -  Ro  ft  r,  -  R0  ft  >/l  +  mSf  (79) 

•=1  i»i 

*,  —  *o  -  ft  M.*.  —  +o  -  ft  Mi  w»‘1/2  tan*l(i,VnT) 

i— 1  i— 1 

The  convergence  of  the  iterations  and  the  efficient  implementation  of  the  rotations  depend 
critically  on  the  predetermined  choice  of  Sf.  Each  type  of  rotation,  (m  =  -1,0,+  1)  has  a 
different  predetermined  set  of  positive  increments  (£<)  which  specify  fixed  radial  and  angular 
increments  (rit  0<)  from  (79).  Within  the  domain  of  convergence  (limited  by  the  total  possible 
rotation)  constraints  were  developed  on  the  sequence  {0,}  such  that  any  angle  0  could  be  rotated 
to  within  0,_t  of  zero  in  p  steps  (WaltherJ.  This  guarantees  that  the  granularity  of  the  calcula¬ 
tion  (the  angular  resolution)  is  0,_j  in  p  steps.  With  the  proper  choice  of  the  set  of  increments 
{$<},  each  successive  iteration  yields  approximately  one  more  bit  of  accuracy  in  the  final  result. 

The  CORDIC  equation,  (78)  is  augmented  by  an  additional  variable,  z,  that  accumulates 
the  angular  component  of  the  rotation. 

h  "  *«  "  ft 

<-i 


(80) 


The  three  parameters  (x,y,z)  are  manipulated  by  successive  application  of  the  CORDIC  equation. 
The  function  to  be  computed  is  obtained  by  forcing  either  y  or  z  to  zero.  Often  the  initial  value 
of  the  other  variable,  z  or  y  is  zero  or  one.  The  sign  of  the  rotation,  p,  is  chosen  at  each  itera¬ 
tion  to  move  the  desired  parameter  towards  zero.  Rotation  operations  are  obtained  by  forcing  zf 
to  zero,  and  vectoring  operations  result  if  y,  is  forced  to  zero. 

The  interpretation  of  (78)  as  a  rotation  on  a  circle,  line  and  hyperbola  can  be  seen  from  a 
graphical  perspective.  For  example,  a  circular  rotation  (m=l)  of  a  vector  \zt  y,| r  into  the  vec¬ 
tor  [z,  0  ]r  computes  lza~lyt/zt.  The  direction  of  rotation  is  chosen  at  each  iteration  to 
force  y,-  closer  to  zero.  The  sequence  of  small  angular  step,  {^},  predetermined  by  {5,},  is  accu¬ 
mulated  with  appropriate  sign  in  zt,  giving  the  answer.  Fig.  8a  indicates  how  this  rotation 
proceeds.  The  radius  of  the  circle  increases  a  predetermined  amount  rt  with  each  rotation  step. 
It  is  not  necessary  to  account  for  this  change  in  radius  when  computing  tan  _1y,/z#.  However  the 
value  of  xf  has  become  a  scaled  square  root  where  the  scale  factor  is  known  in  advance. 

*,  —  y/*S  +  Vo  •  ft  r, 

»= 1 

The  three  input/output  box  of  Fig.  8a  is  used  to  describe  the  function  evaluated  by  the  CORDIC 
rotation.  With  a  nonzero  initial  value  of  z0,  forcing  z,  to  zero  generates  sin(z0)  and  cos(z0),  see 
Fig  8a. 

For  a  rotation  on  a  line,  the  radial  component  is  always  one  and  the  angle  component  is 
interpreted  as  the  y<  value.  The  increments  become  r(  **  l  and  »  6, .  The  result  of  applying 
the  CORDIC  operation  is  shown  in  Fig.  8b.  Multiplication  and  division  can  be  calculated  this 
way. 

The  hyperbolic  rotation  computes  rink,  eo$h,  trctanh  and  square  roots,  see  Fig.  8c.  The 
surface  of  rotation  is  the  set  of  point  that  is  a  constant  distance  vV-y*  from  the  origin.  This 
hyperbola  moves  a  fixed  amount,  closer  to  the  origin  after  each  CORDIC  operation. 


8.2  Lattice  Filtering  by  Rotations 

The  square  root  normalized  lattice  equations  have  a  natural  interpretation  as  rotations.  The 
three  recursive  equations  can  be  efficiently  realized  using  CORDIC  arithmetic  [Ahmed  et  al, 
1981(1),  1981(2)].  The  structure  of  the  SQNLSL  algorithm  suggests  an  implementation  via  rota¬ 
tions.  Since  i>,  p,  and  y  always  have  magnitudes  less  than  one,  they  can  be  interpreted  as  cosines 
of  angles.  Furthermore,  if  x  =»  cot{0,),  then  the  complement  of  x  is  xe  =  \/l  -  x1  =  sin(0s). 
The  SQNLSL  equations  can  be  written  using  a  compact  notation. 

Pi+  i.t  —  Vl-vf.r  v/1  -  li.r-i  Pi+i.r-i  +  Vi.Tli.T- 1  -*■  P+  “  P*neP  +  vri 


vi+  \,T 


Vj,T  -  Pi+  i,r>l«.r-i 
\/l  -  Pi>  i.  r  \/l  -  *7?,  r-i 


*  {v-p+y)l  p'+ye 


Vi+1,T 


Vi.T- 1  ~  Pt+l.Tvi,T 
Vi-Phi.r  V1  ~  vi,T 


1+  =  {l  ~  P+v)  l  P>' 


For  notational  convenience,  the  following  abbreviations  were  used. 

P  ™  Pi+i,r-i  >  P+  “  Pi+i.r 
*  "i.r  .  *V  *  f'i+l.T 
*1  ’li.r-i  »  *7+  *  *7i+i.r 


The  SQNLSL  update  equations  can  be  written  almost  entirely  in  a  single  matrix  equation  (81) 
using  the  rotation  matrices  for  v  and  y. 


v *  v  \p  Ol  I  if*  -y  —  v'y'p  +  vy  vye  -  pyv'  I  —  I  P+  P+v+ 
-v  ve  Lo  lj  y *  v'y-pvy'  v'y'  +  pvy  J  |^+if+  * 


(81) 


On  the  left  hand  side  of  (81),  the  first  matrix  performs  a  rotation  by  0,  *■  cos_,(» /)  and  the  third 
matrix  rotates  by  0f  *■  cos_,(if).  The  result  is  the  complete  update  for  p  and  partial  updates  for  u 
and  q  and  a  term  (*)  of  no  interest.  The  updates  for  v  and  y  are  completed  by  dividing  by  p%. 
This  matrix  equation  (81)  is  directly  realizable  using  CORDIC  operations. 


The  implementation  of  the  SQNLSL  algorithm  in  an  integrated  circuit  proposes  using  two 
processors  in  parallel,  each  executing  sequentially  five  functions,  (Ahmed  et  ai,  1981(1),  1981(2)]. 


The  computation  proceeds  from  (81)  in  three  step;  an  angular  rotation  by  the  multiplication  of 
p  matrix  by  the  r\  rotation  matrix,  and  the  divisions  by  pi .  The  sequence  of  CORDIC  operations 
shown  in  Fig.  10  begins  with  p,  q  and  v  and  computes  and  .  The  functional  elements 

use  the  notation  of  Fig.  9.  The  rotation  angle  0„  is  calculated  as  using  a  circular 

CORDIC  operation  by  processor  2  in  time  slots  1  and  2.  Rotating  the  p  matrix  by  0„  is  computed 
as  two  multiplications  (linear  CORDIC):  p  t)‘  by  processor  1  during  time  slots  1  and  2,  p  t\  by 
processor  2  in  time  slot  3.  These  quantities  are  rotated  by  8„  using  circular  CORDIC  operations 
in  processor  1  (time  slot  3)  and  processor  2  (time  slot  4).  This  generates  the  p  update  and  partial 
updates  for  v  and  rj.  In  time  slot  5,  the  processors  generate  the  updates  for  u  and  q  by  dividing 
the  earlier  results  by  pi .  The  signals  that  flow  farther  than  adjacent  time  slots  must  be  held  in 
temporary  buffers.  Each  CORDIC  operation  uses  16  iterations  and  results  in  almost  16  bits  of 
accuracy.  The  integrated  circuit  could  perform  the  SQNLSL  algorithm  of  tenth  order  on  an  8 
KHi.  sampled  signal  in  real  time.  This  assumes  standard  integrated  circuit  design  rules  to  gen* 
erate  a  moderate  size  chip  running  at  a  20  MHz.  clock  rate. 
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Fig.  8c  CORDIC  Rotation  on  a  Hyperbola 


0.  SIMULATIONS  AND  APPLICATIONS 


SIMULATED  SIGNALS 

Simulated  signals  with  different  characteristics  were  generated  and  the  response  of  the  lat¬ 
tice  estimation  algorithms  noted.  The  simulated  signals  included  white  noise  drive  autoregressive 
processes  that  were  (1)  stationary,  (2)  had  linearly  time  varying  coefficients,  (3)  had  step  changes 
in  coefficients,  and  (4)  had  impulse  excitation  at  the  step  changes  in  coefficients. 

For  a  stationary  autoregressive  process,  the  convergence  of  the  least  squares  lattice  method 
is  shown  in  Fig.  10.  An  eighth  order  fixed  coefficient  lattice  filter  driven  by  white  noise  generated 
the  simulation  data.  The  LSL  algorithm  with  X— .99  was  used  to  compute  the  reflection 
coefficients.  The  first  reflection  coefficient  converged  in  less  than  50  samples  and  the  first  four 
reflection  coefficients  were  near  their  correct  values  after  150  samples.  Higher  order  reflection 
coefficients  approached  their  correct  values  after  250  samples. 

When  the  simulated  data  was  generated  by  a  white  noise  driven  second  order  lattice  with 
linearly  time  varying  coefficients,  the  adaptive  nature  was  be  seen  in  Fig.  11.  The  two  reflection 
coefficient  estimates  followed  the  actual  parameter  values.  However,  there  was  an  increase  in  the 
variance  of  the  estimate  as  the  reflection  coefficients  approached  zero  or  as  the  coefficient  index 
increased.  The  previous  experiment  was  repeated  with  piecewise  constant  coefficients  to  generate 
the  simulation  data.  The  estimated  reflection  coefficient  trajectory  did  not  indicate  that  the 
model  had  step  changes  in  the  coefficients,  see  Fig.  12. 

The  effect  of  the  optimal  weighting  function  7  was  seen  when  the  simulated  data  was  gen* 
crated  by  the  same  lattice  with  step  changes  in  coefficients  but  also  had  a  periodic  impulse  added 
to  the  white  noise  driving  process  at  the  instant  of  coefficient  change.  The  presence  of  the 
impulses  caused  the  estimates  to  readjust  quickly  to  the  new  piecewise  constant  values,  see  Fig. 
13.  The  impulse  caused  a  sudden  increase  in  the  7  which  allowed  the  estimates  to  focus  on  the 
new  signal  characteristics.  Once  the  effect  of  the  impulse  has  passed,  the  7  decreased  so  that  con¬ 
vergence  could  take  place. 


APPLICATION  FOR  SPEECH  ANALYSIS 


The  moat  extensive  use  of  the  lattice  filter  has  been  for  speech  processing  applications, 
including  speech  compression  systems  and  stored  vocabulary  speech  synthesis  chips.  Reflection 
coefficients  and  the  lattice  filter  are  well  suited  for  speech  processing  for  many  reasons;  the  rela¬ 
tionship  to  the  acoustical  tube  model  of  the  vocal  tract,  the  advantageous  quantization  properties 
of  the  reflection  coefficients,  the  finite  word  length  arithmetic  properties  of  the  lattice  filter,  and 
the  slowly  time  varying  nature  of  the  reflection  coefficients  across  speech  sounds  (making  them 
amenable  to  interpolation). 

Linear  Predictive  Coding  (LPC)  is  a  technique  that  has  been  used  widely  for  low  bit  rate 
speech  coding  and  fixed  vocabulary  speech  synthesis.  LPC  uses  a  vocal  tract  model,  the  lattice 
filter  parameterized  by  reflection  coefficients  and  an  excitation  model,  periodic  pulses  for  sounds 
produced  by  vocal  chord  oscillation  (eg.  voweb)  and  white  noise  for  hiss  sounds.  Short  time  seg¬ 
ments  of  speech,  typically  20  milliseconds,  are  characterized  by  eight  to  ten  reflection  coefficients, 
the  pulse  period  (zero  for  noise),  and  an  energy  term.  All  of  the  parameters  can  be  quantized  to  a 
total  of  48  bit  per  20  millisecond  interval.  Using  this  compact  description  of  sounds,  speech  syn¬ 
thesis  integrated  circuit  have  been  developed  that  generate  understandable  speech  using  parame¬ 
ters  store  in  read  only  memory. 

Analyzing  a  spoken  vowel  sound  by  the  LSL  algorithm  shows  the  properties  of  the  likeli¬ 
hood  variable.  The  time  waveform,  Fig.  14a  clearly  shows  the  periodic  nature  of  this  vowel.  The 
LSL  algorithm  applied  to  this  sound  produces  the  forward  prediction  error  shown  in  Fig.  14b. 
This  relatively  stationary  sound  produced  fairly  constant  reflection  coefficients  (after  conver¬ 
gence),  see  Fig.  14c.  The  periodic  jumps  seen  in  all  five  reflection  coefficients  are  due  to  the 
inflaence  of  the  periodic  opening  of  the  vocal  chords.  The  likelihood  variable,  7  usually  b  small 
bat  increases  when  these  openings  occur,  Fig  14d.  When  the  vocal  chords  open,  a  sudden  pulse  of 
air  excites  the  vocal  tract  which  the  likelihood  variable  interprets  as  a  change  in  the  structure  of 
the  signal.  Determining  the  periodicity  of  these  openings,  called  the  pitch  period  is  necessary  of 
the  LPC  speech  model.  Pitch  pubes  can  be  located  directly  from  the  prediction  errors  bat  do  not 


always  give  accurate  results.  The  periodicity  is  evident  but  not  easily  extracted  from  Fig.  14b. 
By  combines  the  derivative  of  the  likelihood  variable  7  with  the  prediction  error  sequence,  a  more 
easily  discernible  spike  is  generated  at  the  onset  of  vocal  chord  oscillation,  see  Fig.  14e.  This 
technique  has  been  proposed  as  a  pitch  estimation  method  [Lee  and  Morf,  1980]. 

Since  the  recursive  exact  least  squares  lattice  algorithms  can  track  quickly  changing  spectral 
characteristics,  they  can  be  used  to  differentiate  the  nature  of  transitional  sounds  [Turner,  1982], 
By  exponential  weighting  of  past  data,  the  current  estimate  reflects  the  short  time  signal  charac¬ 
teristics.  The  beginnings  of  the  words  ’bid’  and  'did'  spoken  by  the  same  male  speaker  are  shown 
in  Fig.  15.  When  analyzed  by  the  SQNLSL  with  X  =  .98,  the  reflection  coefficients  for  the  begin¬ 
ning  of  each  word  follow  different  trajectories  corresponding  to  the  different  consonants.  How¬ 
ever,  during  the  later  vowel  portion,  the  values  are  more  similar,  see  Fig.  15b,c,e,f.  The  transi¬ 
tional  part  of  the  sounds  is  emphasized  but  the  effects  of  the  pitch  pulses  are  also  seen.  The  abil¬ 
ity  of  the  SQNLSL  to  differentiate  these  types  of  sounds  may  be  useful  in  a  phoneme  based 
speech  recognition  systems. 


APPLICATION  FOR  CHANNEL  EQUALIZATION 


The  adaptive  lattice  fitter  offers  substantial  advantages  for  channel  equalization  where  the 
orthogonalizing  properties  and  fast  tracking  characteristics  are  important.  Tapped  delay  line 
adaptive  gradient  equalizers,  although  simple  to  implement,  have  a  rate  of  convergence  that 
depends  on  the  ratio  of  the  largest  to  smallest  eigenvalues  of  the  channel  correlation  matrix  [Git- 
lin  et  a),  1973],  Self-orthogonalizing  techniques  have  been  proposed  by  [Gitlin  and  Magee,  1977] 
and  in  lattice  form  by  (Griffiths,  1977,  Griffiths  and  Medaugh,  1979(2)].  The  gradient  lattice 
equalizer  [Satorius  and  Alexander,  1979(2)]  and  the  LSL  equalizer  (see  Section  6)  [Satorius  and 
Pack,  1981]  have  been  shown  to  provide  very  fast  convergence.  The  lattice  filter  equalizers 
demonstrated  fast  convergence  that  was  independent  of  the  channels  eigenvalues  disparity  ratios, 
see  Fig.  16  from  (Satorius  and  Pack,  1981],  Two  simulated  data  channels  with  correlation 
matrices  of  eigenvalues  disparity  ratios  (ratio  of  largest  to  smallest  eigenvalues)  of  11  and  21 
respectively  were  studied.  An  11  tap  equalizer  was  implemented  using  the  LMS  gradient  algo¬ 
rithm,  the  gradient  lattice  algorithm  and  the  LSL  algorithm.  The  gradient  tapped  delay  line 
equaliser  had  considerably  slower  convergence  that  depended  on  the  eigenvalue  ratio.  The  LSL 
equaliser  converged  in  both  cases  in  approximately  40  iterations  while  the  adaptive  lattice  equal¬ 
iser  required  approximately  120  iterations  to  converge. 


Figure  16 


APPLICATION  FOR  ELECTROENCEPHALOGRAPHIC  ANALYSIS 


Electroencephalographic  (EEG)  data  analyzed  by  autoregressive  modeling  can  provides  a 
better  summary  of  EEG  spectral  information  than  frequency  domain  techniques,  such  as  FFTs. 
The  reflection  coefficients  from  the  SQNLSL  algorithm  were  studies  to  detect  subtle  changes  in 
brain  states  as  observed  in  EEG  activity  (Redington  and  Turner].  The  data  obtained  from  the 
left  central  EEG  (Cl)  response  of  an  adult  human  subject  monitored  during  sleep  onset  (sampled 
60  times  a  second)  is  shown  in  Fig.  17a.  A  large  change  in  activity  appears  near  the  beginning  of 
the  raw  EEG  data  trace  and  is  apparent  in  the  reflection  coefficients,  Fig.  17b,c,d.  A  second 
change  in  activity  near  the  end  of  the  trace  is  barely  noticeable  in  the  raw  data;  yet,  it  is  easily 
recognized  in  the  activity  of  the  higher  order  reflection  coefficient.  The  changes  in  reflection 
coefficients  may  reflect  physiological  transitions  and  provide  a  means  of  inferring  presence  or 
sequence  of  EEG  brain  states. 


10.  COMMENTS  AND  CONCLUSIONS 


The  lattice  filter  structure  as  a  realization  of  a  digital  transfer  function  has  several  advan¬ 
tages;  it  is  a  cascade  of  identical  sections,  has  a  general  insensitivity  to  round-off  noise,  and  the 
reflection  coefficients  can  be  related  to  physical  processes.  The  physical  interpretation  of 
reflection  coefficients  gives  them  intuitive  appeal,  particularly  for  speech  signals.  For  adaptive 
estimation,  the  lattice  structure  is  the  natural  form  for  an  efficient  solution  to  recursive  least 
squares  problems.  Lattice  filters  provide  an  orthogonalization  or  decoupling  of  the  states  of  the 
input  process.  The  stability  of  an  all  pole  model  when  expressed  in  lattice  form  can  be  deter¬ 
mined  by  inspection. 

The  real  advantage  of  the  lattice  structure  lies  in  adaptive  estimation  and  filtering.  An  N 
stage  lattice  filter  automatically  generates  all  the  outputs  which  would  be  generated  by  N 
different  TDL  filters  with  lengths  from  1  to  N.  This  allows  dynamic  assignment  of  any  filter 
length  which  proves  most  effective  at  any  instant  of  adaptive  processing.  When  compared  to  the 
simpler  adaptive  transversal  filter,  the  lattice  filter  has  superior  convergence  properties  and 
reduced  sensitivity  to  finite  wordlength  effects. 

Recursive  lattice  estimation  algorithms  allow  the  exact  least  squares  solution  to  be 
efficiently  updated  for  each  new  time  sample.  The  structure  of  this  exact  recursive  approach  is 
similar  to  the  gradient  lattice  techniques;  however  an  optimal  gain  is  calculated  at  every  time 
sample.  This  optimal  recursive  solution  has  a  complexity  that  is  only  slightly  more  that  the  gra¬ 
dient  lattice  solution.  Consequently,  the  LSL  and  SQNLSL  algorithms  achieve  extremely  fast  ini¬ 
tial  convergence  and  can  track  quickly  time  varying  parameters.  The  SQNLSL  has  a  very  com¬ 
pact  notation  and  normalizes  all  signals  to  unit  variance  at  each  stage.  A  single  integrated  circuit 
to  execute  this  algorithm  has  been  proposed. 

However,  as  with  all  adaptive  estimation  procedures,  there  are  various  trade-offs  to  be  made. 
The  lattice  structure  involves  more  computation  and  is  conceptually  more  complicated  than  the 
tapped  delay  line  structure  but  has  better  convergence  properties.  The  recursive  least  squares  I  at- 


tice  offers  even  better  convergence  than  the  gradient  lattice,  but  again  it  is  slightly  more  complex. 
For  example,  in  the  definition  of  the  two  reflection  coefficients,  there  is  a  difference  in  the  time 
subscripts  of  the  normalising  covariances.  In  the  stationary  case,  these  terms  are  identical  but  in 
the  LSL  the  difference  is  critical;  in  general  the  algorithm  will  fail  if  this  difference  is  overlooked 
[Satorius  and  Shensa,  1980).  The  SQNLSL  allows  a  very  short  time  constant  to  be  applied  to  the 
sampled  data  so  that  the  quickly  time  varying  nature  of  the  signal  can  be  tracked.  However, 
attempting  to  track  transient  speech  sounds  also  tracks  the  pitch  excitation  signal  which  was  not 
of  interest.  For  processes  that  tend  toward  stationarity,  the  convergence  properties  of  the  gra¬ 
dient  lattice  and  LSL  lattice  are  similar  (Honjg,  1983]. 

Many  extensions  to  the  basic  recursive  least  squares  algorithm  have  been  developed. 
Reviews  of  least  squares  adaptive  lattice  filtering  can  be  found  in  [Satorius  and  Shensa,  1980, 
Friedlander,  1982(3)).  Recursive  ladder  algorithms  for  ARMA  modeling  have  been  presented  in 
[Lee  et  al,  1982).  The  SQNLSL  algorithm  has  been  extended  from  the  "pre-windowed”  data  case 
presented  here  to  the  "Covariance”  data  case  in  [Porat  et  al,  1982).  The  problem  of  system 
identification  has  been  addressed  in  [Porat  and  Kailath,  1983|.  A  review  of  lattice  filters  for  nons¬ 
tationary  processes  was  presented  in  [Kailath,  1982). 

There  are  other  means  to  implement  the  lattice  filter  structure  for  estimation.  The  order- 
update  recursions  can  also  be  obtained  by  using  a  Cholesky  decomposition  of  the  covariance 
matrix  [Dickinson,  1979(1),  Dickinson  and  Turner,  1979(2),  Klein  and  Dickinson).  Alternatively, 
since  a  reflection  coefficient  is  similar  to  a  correlation  coefficient,  computationally  simple  tech¬ 
niques  to  estimate  correlation  coefficients  can  be  applied  to  determining  the  reflection  coefficients. 
Since  the  correlation  of  Gaussian  random  variables  is  related  to  the  correlation  of  the  hardlimited 
variables  by  an  ARCSIN  relationship,  a  very  simple  reflection  coefficient  approximation  technique 
is  possible  (Turner  et  al,  1980).  This  algorithm  requires  only  a  count  of  polarity  changes  in  the 
prediction  errors  to  estimate  the  reflection  coefficients  (assuming  zero  mean  unit  variance  Gaus¬ 
sian  signals). 

Overall  the  adaptive  lattice  filter  offers  a  compact  algorithm  for  obtaining  quickly  converg- 


ing  estimates.  The  properties  of  the  lattice  filter  and  reflection  coefficients  motivate  their  use  in 
many  practical  situations. 
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